J X**' 


AD-A009  716 


ANALYSIS  OF  FLOW  AND  HEAT  TRANSFER  IN  SMALL  ARMS 
Edward  V.  McAssey,  Jr. 

Villanova  Uni  'ersity 


Prepared  for: 

Army  Research  Office 
Frankford  Arsenal 


1 April  1975 


DISTRIBUTED  BY: 


IIaKaimI  i— ^ m la—  alii* 

MUVIWB  IGWHIIIMI  111  I III  MMUIIII  OVIVIWV 

U.  S.  DEPARTMENT  OF  COMMERCE 


m**Ksmctirt9TQ* 


iU'  v'^rvnntt^rSk 


Unclassified 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  D*ta  Entered) 

REPORT  DOCUMENTATION  PAGE 


I REPORT  NUMBER 


2.  GOVT  ACCESSION  NO. 


I?EAD  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 
3.  RECIPIENT'S  CATALOG  NUMBER 


Project  II 286 


4.  TITLE  (and  Subtitle) 

Analysis  of  Flow  and  Heat  Transfer  in  Small 
Arms 


ID-  A 00  9 


5.  TYPE  OF  REPORT  & PERIOD  COVERED 

Final 

1 Nov.  73  to  1 April  75 

6,  PERFORMING  ORG.  REPORT  NUMBER 


I 7.  AUTHOR^; 


6.  CONTRACT  OR  GRANT  NUMBERfa) 


Edward  V.  McAssey,  Jr. 


DAHC04-74-G-0044 


9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Mechanical  Engineering  Department 
Villanova,  University,  Villanova,  Pa.  19085 

II.  CONTROLLING  OFFICE  N AME  AN  O ADDRESS  ~~  ~™ 

U.  S.  Army  Research  Office 
Box  CM,  Duke  Station 
Durham,  North  Carolina  27706 

14.  MONITORING  ‘.GENCY  NAME  ft  ADORESSfJf  dlllerent  Iron  Controlling  Otllce) 

Frankford  Arsenal 
Bridge  & Tacony  St. 

Philadelphia,  Pa.  19137 

16.  DISTRIBUTION  STATEMENT  (ol  thla  Report) 


to"  PROGRAM  ELEMENT.  PPOJECT,  TASK 
AREA  ft  WORK  UNIT  NUMBERS 


P-12080-RTL 


12.  REPORT  DATE 

1 Aprjj.  1975 

13.  n/mbeKof  pagi 


mbeAof  pages 

£ 5" 


15.  SECURITY  CLASS,  (ol  title  report) 

Unclassified 

“16*.’  DECLASSIFICATION' DOWNGRADING 
SCHEDULE  jj* 


Approved  for  public  release  distribution  unlimited. 


17.  DISTRIBUTION  STATEMENT  (ol  the  abatraet  entered  In  Block  20,  II  dlllerent  Irom  Report) 


!G.  SUPPLEMENTARY  NOTES 


ID.  KEY  WOROS  (Continue  on  teverae  side  if  necessary  and  Identify  by  Vocfr  number) 

Transient  Flow,  Heat  Transfer,  Flow  in  Small  Arms,  Flame  Front 


A05TRACT  (Contlnuo  oil  ravarae  aide  If  nacoaaety  and  Identify  by  block  number) 

A computer  model  of  the  pressure  distribution  and  flow  characteristics  in 
<5mal  ] arms  barrels  has  been  developed.  This  mode'.  «hows  the  senai  tHv  i tv  of 
the  pressure  history  and  launch  velocity  of  surface  concentration  density  of 
deterrent,  flame  front  velocity,  and  projectile  base  heating. 


Reproduc'd  by 

NATIONAL  TECHNICAL 
INFORMATION  SERVICE 

US  Department  ol  Commerce 
Springfield,  VA.  22151 


PRICES  SlMJECT  TO  CHANGE 


I FAfT! 73  1473  COITION  OF  1 NOV  65  IS  OBSOLETE 


Unclassified 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Drte  Entered) 


rjf-ar?.’?  VJ>  v~*ri 


I 0 


1 If 


fj  S: 

ft  at 

*- 

&: 


I 


I 


r 


ANALYSIS  OF  FLOW  AND  HEAT  TRANSFER 
IN  SMALL  ARMS 


Edward  V.  McAssey,  Jr. 
April  1,  1975 


Sponsored  by  U.S.  Army  Research  Office 
Durham*  N.  0. » Grunt  DAHCQ4-74-G-0044 
Project  No.  286 


Mechanical  Engineering  Department 
Villanova  University 


Approved  for  Public  Release; 
Distribution  Unlimited 


H 


Miiil  f f-n  to  rf» » i afeaa Sills  a aae-  ■■ » *■  ? s.  .-r 


I 


Abstract 

A computer  model  of  the  pressure  distribution  and  flow  characteristics  in 
small  arms  barrels  has  been  developed.  This  model  shows  the  sensitivity 
of  the  pressure  history  and  launch  velocity  to  surface  concentration  density 
of  deterrent,  flame  front  velocity,  and  projectile  base  heating. 


ALL  FINDINGS  IN  THIS  REPORT  ARE  NOT  TO  B£  CONSTRUED  AS  AN  OFFICIAL 
DEPARTMENT  OF  THE  ARMY  POSITION,  UNLESS  SO  DESIGNATED  BY  OTHER  AUTHORIZED 
DOCUMENTS. 
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NOMENCLATURE 

A cross-sectional  area  of  gun  boro,  ft^ 

a speed  of  sound,  ft /sec 

Cp  specific  heat  of  propellant  gas  at  constant  pressure,  BTU/lbm-°R 

c.y  specific  heat  of  propellant  j as  at  constant  volume,  BTU/lbm-°R 

0 diameter  of  gun  bore,  ft 

D0  initial  diameter  of  propellant  grain  (weighted  average),  ft 

Egas  energy  per  unit  mass  of  gas,  BTU/slug 

Ej  drag  force  on  propellent  particles  per  unit  volume,  lbf/ft^ 

g acceleration  of  gravity,  ft/sec^ 

H enthalpy,  BTU/slug 

J mechanical  equivalent  of  heat,  ft-lbf/BTU 

1 distance  burned  normal  to  propellant  grain  surface,  ft 

mg  rate  of  mass  of  gas  generated  per  unit  length,  slug/ft-sec 

p gas  pressure,  lbf/ft^ 

q heat  loss  through  barrel  wall,  BTU/ft  -sec 

R gas  constant  for  propellant  gas,  f t-lbf /slug-°R 

R disk  radius  fo  an  individual  grain,  ft 

Ri  one  half  of  web  thickness  of  individual  grain,  ft 

r linear  ourning  rate  of  propellant,  ft/sec 

s specific  entropy,  BTU/lbm-°K 

Sp  surface  area  of  an  individual  grain,  ft^ 

Tg  gas  temperature,  °R 

t time,  sec 

u velocity  of  gas,  f ' /sec. 

b 

u.  Vvtui  icy  of  1 luuL  ^uiliclcb*  fi/s»cc 


% 


volume  of  an  individual  grain,  ft2 

powder  potential,  cvT(i),  BTU/lbm 

weight  of  projectile,  lbm 

particle  volume  fraction,  dimensionless 

number  of  particle  grains  per  unit  volume,  ft--* 

gas  density,  slug/ft^ 

propellant  density,  slug/ft-* 

wall  shearing  stress,  lbf/ft2  or  slug/ft-sec2 
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Introduction 

Investigations  into  propellant  performance,  barrel  materials  and 
weapon  kinematics  require  accurate  description  of  the  weapon  interior 
ballistics.  The  present  study  was  undertaken  to  determine  the  effect  of 
heat  transfer  on  interior  ballistics.  In  examining  the  interior  ballistics 
of  small  arms,  the  effect  of  a number  of  parameters  must  be  included.  These 
parameters  include  chamber  geometry,  pr^a1 lant  geometry,  propellent  com- 
position and  burning  rate. 

The  heat  transfer  rates  in  small  arms  are  very  high  and  are  considered 
to  be  a major  cause  of  barrel  errosion.  Accurate  description  of  these  heat 
transfer  rates  is  necessary  for  analytical  predication  of  pressure  and 
velocity  histories  in  small  arms. 

The  present  study  uses  the  method  of  characteristics  to  solve  the 
governing  equations.  The  model  assumes  the  problem  to  be  one  dimensional 
and  transient.  Assuming  a fixed  propellent  description,  the  effect  of  var- 
ious heat  transfer  models  has  been  examined.  The  performance  standard  is 
test  data  from  reference  1 presented  in  figure  1. 
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literature  Survey 

Reference  (2)  presents  a ballistic  model  developed  from  earlier  work 
which  was  modified  to  include  the  effect  of  non-uniform  propellant  compo- 
sition and  non-simul taneous  ignition.  In  this  reference  the  propellant 
is  represented  by  a composite  of  9 different  composition  propellants. 

This  variable  propellant  array  allows  the  handling  of  different  siie  grains 
and  different  deterrent  concentration  densities.  Sample  results  are  in- 
cluded in  reference  (2)  which  shows  good  agreement  with  test. 

Reference  (3)  presents  an  analytical  investigation  of  the  flow  in  gun 
barrels  with  particular  emphasis  on  the  hear,  transfer.  In  addition  to  the 
transient  nature  of  the  flow,  reference  (3)  considers  the  problem  posed  by 
the  lack  of  a boundary  layer  at  the  base  of  the  projectile.  This  condition 
occurs  because  in  the  projectile  based  region  the  gas  velocity  is  equal  to 
the  projectile  velocity.  This  condition  means  that  the  Reynolds  analogy 
commonly  used  in  heat  transfer  cannot  be  applied  to  this  case.  However, 
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the  reference  does  propose  an  expression  for  heat  transfer  coefficient  in 
the  projectile  base  regioi.  The  resulting  convective  coefficients  are  be- 
tween 5 and  20  Btu/sec  f t^  °F.  Comparisons  between  experiement  and  pre- 
diction presented  in  reference  (3)  indicate  that  even  these  high  coefficients 
underestimate  the  actual  results.  However,  the  results  of  reference  (3) 
represent  the  best  available  procedure  for  calculating  the  heating  at  the 
base  of  the  projectile. 

Reference  (4)  presented  an  analytical  study  of  interior  ballistics  in 
which  the  governing  equations  were  solved  using  a finite  difference  procedure. 
The  geometric  dimensions  used  in  reference  (4)  did  not  correspond  to  standard 
small  arms.  In  addition,  the  analysis  was  limited  to  stationary  propellant 
and  to  propellant  moving  at  the  same  velocity  as  the  gas.  In  the  area  of 
heat  transfer  reference  (4)  also  obtained  very  high  heat  transfer  coefficients 
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10  3tu/sec  fc2  °F) . It  is  worthy  of  note  that  for  cases  with  and 
without  heat  transfer  and  wall  friction  the  ballistic  performance  was  not 
changed  significantly.  The  present  study  has  found  that  heat  transfer  par- 
ticularly at  the  projectile  base,  is  very  significant. 

Reference  (5)  presents  an  extensive  study  of  the  flame  front  movement 
through  a porous  bed  propellant.  The  analysis  is  presented  for  a venting 
bomb  type  device.  However,  the  results  for  flame  propagation  can  be  applied 
to  the  small  arms  case.  Indeed,  the  functional  relationship  between  flame 
velocity  and  pressure  can  be  obtained  from  the  data  presented  in  reference  (5). 
Reference  (6)  considers  the  same  problem  as  reference  (5)  but  applied  to  the 
small  arms  case.  The  results  of  reference  (6)  have  been  used  in  the  present 
study  to  determine  the  progress  of  the  flame  front  through  the  propellant  bed. 
Reference  (7)  presents  an  iterative  model  for  the  interior  ballistic  problem 
in  which  parameters  affecting  the  burning  rate  were  changed  until  the  calcu- 
lated results  agreed  with  test  data.  The  caliber  used  in  reference  (7)  is 
the  same  as  the  present  study,  (i.e.  5.56  MM);  however,  the  experimental  and 
analytical  results  are  considerably  different.  This  is  most  likely  due  to 
different  charge  weight  and  propellant  type.  In  addition,  a number  of  the  loss 
factors  considered  in  the  present  study  are  not  included. 

Reference  (8)  is  the  analytical  basis  for  the  present  study.  This 
reference  presents  the  model  which  has  been  used  and  also  outlines  the  cal- 
culation procedure  for  the  method  of  characteristics.  The  nomenclature  of 
reference  (8)  has  been  adopted  in  the  present  report.  The  present  investiga- 

* 

tion  has  expended  considerable  effort  in  programming  reference  (8).  Refer- 
ence (8)  describes  the  interior  ballistic  problem  as  a one  dimensional 
transient  flow  problem.  In  the  model  the  propellant  is  assumed  to  ignite 
instantaneously  and  burn  at  a rate  proportional  to  the  gas  pressure.  The 
geometry  of  the  propellant  is  uniform  and  based  upon  reference  (10). 
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Although  heat  transfer  between  the  gas  and  propellant  is  neglected,  drag 
between  these  two  is  considered.  The  transformation  of  the  governing 
equations  from  the  physical  domain  to  the  characteristic  domain  is  described 
in  detail  in  reference  (9).  The  heat  transfer  to  the  wall  and  the  wall 
friction  are  included  in  the  model  but  the  detailed  descriptions  are 
developed  in  the  present  study.  In  addition,  the  effect  of  chamberage 
and  projectile  resistive  forces  are  also  included  in  the  model  of  reference 
(9). 
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Analysis 

T',  determine  the  effect  of  heat  transfer  and  to  evaluate  potential 
heat  transfer  relationships,  an  accurate  model  of  the  flow  field  must  be 
developed.  Considerable  effort  during  the  present  study  was  devoted  to 
developing  a model  which  could  reproduce  test  data.  In  addition,  during 
the  course  of  this  model  development  it  became  apparent  that  the  assumption 
of  uniform  ignition  was  not  adequate  and  therefore  a more  realistic  flame, 
front  model  had  to  be  developed.  This  section  describes  the  flow  field 
flame  front  analysis,  and  heat  transfer  models. 

The  flow  is  two  phase,  non-steady,  including  heat  and  mass  addition. 

The  basic  model  describing  the  problem  has  been  developed  in  reference  (8). 
Because  of  the  geometry  present,  the  problem  is  quasi-one  dimensional.  The 
resulting  conservation  equations  form  a set  of  quasi-1 inear  partial  differ- 
ential equations  of  the  first  order  and  hyperbolic  type.  For  this  situation 
the  method  of  characteristics  provides  the  simplest  means  of  solution. 
Reference  (8)  provides  the  details  of  the  transformation  of  the  governing 
equations  to  the  characteristic  form.  In  addition,  reference  (9)  provides 
who  theoretical  foundation  for  the  applications  of  the  method  of  characteris- 
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tic  to  the  present  problem. 


Flow  Field  Governing  Equations 

The  governing  equations  described  in  these  sections  are  based  upon 
the  following  assumptions: 

1.  Flow  is  quasi-one  dimension.il 

2.  The  gas  is  assumed  to  be  ideal 

3.  The  specific  heat  and  molecular  weight  are  dependent  upon 
the  amount  of  deterred  and  non-deterred  propellant  burned 
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4.  All  propellants  grains  are  initially  of  uniform  surface  area,  volume, 
and  radius.  Weighted  mean  values  based  upon  reference  (10)  and  (11) 
are  used. 

5.  Interactions  between  individual  propellant  particles  and  between 
particles  and  the  boundary  are  negligible 

6.  The  thermal  motion  of  the  propellant  particles  does  not  contribute 
to  the  pressure 

7.  The  average  gas  temperature  is  defined  by  the  heat  content  of  the 
gas  phase 

K.  The  propellant  particles  are  spheres  for  drag  calculations 

9.  Heat  dissipation  due  to  viscosity  is  negligible 

10.  Radiation  is  negligible 

11.  Gravity  is  negligible 
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Overall  Gas  Continuity  Equation 
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where  rag  is  the  gas  generated  per  unit  length 
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Propellant  Particle  Continuity  Equ.it  i on 
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where  Tlw  is  shear  stress  at  the  wall  and  Fd  is  particle  drag  per  unit  volume 


Propellant  Particle  Momentum  Equation 
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Overall  Gas  Energy  Equation 
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where  W is  the  heat  content  of  the  propellant  per  pound 
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Overall  Gas  Entropy  Equation 
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In  the  preceding  equations,  the  dependent  gas  variables  are  entropy, 
gas  temperature,  pressure,  gas  density  and  gas  velocity. 

The  next  step  in  the  solution  procedure  is  to  apply  the  method  of 
characteristics.  In  this  approach  the  dependent  variables  become  the 
P and  Q Riemann  variables  where: 

~p=  Z CL  + Uo 

*»«l  * 

Q a -2l.  CL  * Uq 

* 

Solutions  are  determined  along  these  P and  Q characteristics. 

The  characteristics  represent  lines  along  which  infinitesimal  disturbances 
are  pripagated  at  the  local  speed  of  sound  in  the  gas  phase.  Along  each 
characteristic  the  following  relationships  between  X,  t are  true: 
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Equation  (1)  ..nd  (3)  can  be  combined  into  two  equations  for  P and  Q along 
the  respective  character ist ic  using  the  ideal  gas  relations  (for  details  see 
reference  (3))j  qX  _ *2. 

JLPa  H JL»Ol  - VLS 
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yielding  equations  (7>  and  (8). 
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Equation  (6),  (7)  and  (8)  along  with  the  appropriate  ideal  gas  relations 
are  the  governing  equations  for  the  gas  phase.  Equations  (2)  and  (A)  describe 
the  propellant  particle  phase. 

Initial  and  Boundary  Condtions 

At  t“0,  the  first  Q wave  leaves  the  base  of  the  projectile  and  travels 
toward  the  breech.  During  the  time  prior  to  the  arrival  of  the  Q wave  at 
the  breech,  it  is  assumed  that  the  gas  in  the  chamber  is  at  rest  and  that 
the  temperature  is  constant  at  the  flame  temperature,  Tq.  The  pressure  in 
the  region  between  the  projectile  and  breech  will  rise  as  the  Q wave  advances 
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due  to  the  burning  of  propellant.  The  slope  of  this  first  Q wave  is  given 


by  equation  (9) 
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The  Q wave  is  divided  into  a specified  number  of  equally  spaced  nodes. 

The  nodes  are  the  origin  for  the  P wave.  At  each  node  the  gas  temperature 
and  speed  of  sound  are  known,  the  other  parameters  are  given  by  the  following 
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Strict;  the  speed  of  sound  is  constant  along  this  wave,  Q is  also  constant. 
This  procedure  is  different  from  that  described  in  reference  8,  but  it 
yields  the  same  results. 

The  spatial  boundaries  of  the  problem  are  the  breech  (x*0)  and  the  base 
of  the  projectile.  The  breech  is  at  a fixed  condition  with  the  gas  and 
propel!  lnt  velocity  equal  to  zero.  The  location  of  the  projectile  base 
is  determined  by  the  intersection  of  a P wave  and  the  equation  of  motion 
of  the  projectile.  The  two  equations  governing  the  position  of  the 
projectile  are: 
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where  is  the  engraving  force  and  Wt  is  the  weight  of  the  projectile. 

At  the  projectile  base  the  gas  velocity,  propellant  velocity  and  projectile 
velocity  are  all  equal. 


There  are  three  types  of  grid  points,  interior  points,  (a,b,c), 
breech  points  (e,f.g)  and  projectile  base  points  (h,  l,  j).  The  govern- 
ing equations  and  solution  procedures  are  different  for  each  of  these  types 
of  grid  points  and  therefore,  each  will  be  handled  separately. 


I 

I 

I 


I 


w.\ 


p 

u 


fi 

u 


a 


Interior  Points 


Figure  2 shows  a typical  interior  unknown  grid  point  (6  . it  is  on 
a P wave  from  point  (4)  and  a Q wave  from  point  (3) . 
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Figure  2:  Diagram  of  Points  Used  to  Calculate 

Interim  Grid  Points 


Point  (7)  represents  the  location  of  a gas  particle  which  intersects 
point  (6)  and  . he  line  connect!  ig  points  (3)  and  (4) . Point  (5)  is  similar 
to  point  (7)  except  it  represents  a propellant  particle.  The  slope  of  the 
P and  Q characteristics  are  given  by  equations  (12)  and  (13) . 
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Solving  for  t$  and  X& 
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Since  the  values  at  point  6 are  initially  unknown,  an  assumed  value 
of  a given  parameter  is  used  and  then  the  procedure  is  repeated  until  the 
assumed  value  and  calculated  value  agree  to  within  a specified  tolerance. 
The  initial  ; , caption  is  an  arithmetic  average  of  the  v. lues  at  points  (3) 
and  (4). 

Knowing  the  location  of  point  6,  it  is  possible  using  the  procedure 
described  in  reference  (8)  pages  23  through  30  to  generate  the  data 
required  for  the  solution  to  proceed.  For  brevity,  the  details  will  be 
omitted;  however,  certain  key  points  should  be  mentioned.  The  average 
burning  rate  between  points  (5)  and  (6)  given  by: 

S fc  ( l + -V  ' 'VUf'vf  So*^ 


The  distance  burned  normal  to  a particle  surface  at  (6)  is: 
*6  * 15  + r5-6  (tg-t5> 


Knowing  V6  and  «. 5 the  average  particle  volume  and  suiface  area  between 
(5)  and  (6)  can  be  determined 


Vp5,6  = <Vp5<*.5)  + Vp6U6))/2 


Sp5,6  ~ + 8p^(£6})/2 

The  gas  per  unit  length  generated  along  the  propellant  path  (5)  ar.d  (6)  is 
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The  rate  of  change  in  propellant  particle  velocity  between  (5)  and  (6)  is 
given  by  equation  (16) 
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The  rate  of  change  of  propellant  particle  volume  fraction  is  given  by 
equation  (17) 

*2i  « ~ *»>  "5^?  - 2 €.t,i  Uyg,i  ( &4-  &») 

- (17) 

M* 

The  rate  of  entropy  change  for  the  gas  can  now  be  computed  using 
equation  (6).  The  change  in  P and  Q along  the  appropriate  characteristic 
is  given  by  equations  (18)  and  (19) 
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Using  the  calculated  values  for  P^  and  Qg,  new  values  for  the  gas 
velocity  ana  speed  of  sound  can  be  computed 


V * (P6-Q6)/2 
a6  **  <i=I)  (p6  + q6> 


The  entropy  at  point  6 can  be  calculated  from  equation  (24) 
S6  * s?  + (~)  (t6  - t7) 


The  new  values  of  Ugg,  a^,  and  S5  can  be  used  to  calculate  the 
remaining  thermodynamic  properties  based  upon  the  equations  on  page 
34  reference  8. 

In  addition  Ug^  and  are  available  to  obtain  a new  estimation  of 
the  location  of  point  (6)  (i.e.,  t6,  X6).  The  procedure  described 
above  is  then  repeated  until  the  new  values  of  a<j  and  Ug$  agrees  with 
Che  previously  calculated  ones  to  within  a specified  tolerance. 
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Projectile  Base  Point 

The  solution  at  the  projectile  base  is  determined  by  the  intersection 
of  a P wave  and  the  path  of  the  projectile.  Figure  3 presents  a diagram 
of  this  intersection.  The  location  and  thermodynamic  properties  at 
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Figure  3:  Diagram  of  Points  Used  to  Calculate  a 
Projectile  Base  Point 


points  1 and  2 are  known.  This  section  describes  the  procedure  to  find 
point  4.  The  slope  of  the  P wave  from  2 to  4 is  given  by  equation  (25) 


*4" 

*4-  t. 


UVu  + Si 


In  addition  integration  of  the  equation  of  motion  of  the  projectile 
yields  equation  (26) 
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Simultaneous  solution  of  equations  (25)  and  (26)  yields  equation  (27) 

^*4  * " ^ ^ S*-  AAC 
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where 
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Knowing  t4»  X4  can  be  obtained  from  equation  (25).  At  the  base  of  the 
projectile  the  projectile  velocity  equals  the  gas  velocity  (eq.  (28)). 

UV.UV* 


(28) 


Having  established  a first  approximation  for  the  location  and  the  velocity 
at  point  4,  the  change  in  P along  the  characteristic  2-4  can  be  determined. 
The  details  regarding  this  procedure  are  given  in  reference  8,  pages  36 
to  40.  The  rates  of  change  of  P ( S*V/  Si  ) and  P4  are  given  by 


equations  (29)  and  (30) 
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Knowing  P4  and  Ug^,  a new  velocity  of  sound,  34,  can  be  computed  from 
equation  (31) 

flt4  - ( ?4-  (31) 


Using  34  and  S4,  revised  values  for  the  other  thermodynamic  variable  can 
be  computed.  The  procedure  is  then  repeated  until  the  values  of  Ug4  and 
34  for  two  consecutive  calculations  agree  within  a specified  tolerance. 
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After  the  projectile  leaves  the  barrel  the  procedure  is  the  same  as 
described  on  pages  49  and  50  of  reference  8.  It  is  assumed  that  the  gas 
velocity  is  sonic  at  muzzle  exit  and  that  Q waves  after  the  projectile 
leaves,  originate  from  the  muzzle  exit.  Under  these  conditions  X24  is 
fixed  and  the  following  relations  exist: 
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Breech  Points 


The  breech  is  located  at  X»0.  Figure  4 diagrams  the  intersection 
of  a Q wave  with  the  breech. 


Figure  4:  Diagram  of  Points  Used  to  Calculate 


Breech  Point 


The  time  t3  is  given  by  equation  (32) 

t3  = -X2/Ug2  - ^2’  3)  + fc2  ^ 

In  equation  (32)  the  initial  estimate  if  a^  is  made  by  assuming  the 
same  values  as  point  (1).  The  change  in  Q from  2 to  3 is  then  computed 
and  a revised  a3  can  be  obtained  from  equation  (33) 
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The  rate  of  change  of  propellant  volume  fraction  and  of  gas  entropy  are 
given  by  equations  (34)  and  (35) 
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where 
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The  rate  of  change  in  Q along  the  characteristics  from  2 to  3 is  given 
by  equation  (36) 
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The  value  of  Q and  S at  point  3 is  given  by  equation  (37)  and  (30) 
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Using  equation  (33)  and  (37),  a new  value  of  83  can  be  calculated.  In 
addition,  aj  and  S3  can  be  used  to  calculate  the  remaining  thermodynamic 
variables.  The  procedure  is  repeated  until  the  value  at  frot  two 
successive  calculations  agrees  within  a specified  tolerance. 
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The  preceding  procedure  has  been  programmed  and  the  details  of  the  program 
are  presented  in  Appendix  A. 

Flame  Front  Governing  Equations 

Initially  it  was  assumed  that  there  was  instantaneous  ignition  of  the 
entire  propellant  bed.  This  resulted  in  exceseive  breech  pressure  and  low 
projectile  base  pressure.  A better  approximation  is  obtained  by  using  a flame 
front  with  a finite  velocity.  Based  upon  the  results  from  references  (4)  and 
(12),  it  was  determined  that  the  flame  front  velocity  was  proportional  to  the 
pressure  behind  the  front.  Using  data  from  these  references  the  following 
relationship  was  developed 


Vp  » op 
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(39) 


where 


a * 31.5  in/sec/(spi) 
8 * <48 
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The  effect  of  the  flame  front  progressing  through  the  propellant  bed 
can  then  be  determined  by  integrating  from  the  breech  to  the  projectile  base. 
The  pressure  and  temperature  behind  the  flame  front  was  assumed  uniform.  The 
gas  velocity  was  also  assumed  to  be  zero.  For  a given  time  increment  the 
distance  the  flame  front  moves  is  given  by  equation  (40). 
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VF  dt 


(40) 


In  equation  (40)  the  velocity  is  based  upon  the  initial  pressure(at  tf) 
duu  au  uiucu  miui  picobiui:  yui  . iut!  pru^uuure  ib  rcpuatea  until  tnu 
assumed  final  pressure  and  calculated  final  pressure  agree  to  within  5Z.  The 
distance  burned  normal  to  the  propellant  surface  in  the  same  time  period  is 


given  by  equation  (41) 


d XL  - (B  pn  ) (t2-tL) 


(41) 


Using  the  results  of  reference  (10)  the  new  volume  a propellent  grain  resulting 
from  the  surface  recession  given  in  equation  (41)  is 

A,)  * JLl, >-i  +1*1 
RL  » - A,} 

: tn  tmMtRlA  ♦ n*gl*(  RR  + +RO  <42> 

^ *1 


where 

RR  = R0  - Ri 

R0  « equivalent  mean  radius  of  propellant  grain 
R^  * one  half  web 

The  total  gas  generated  in  the  cell  defined  by  dL  from  the  time 
t * 0 to  t s t2  is  given  by  equation  (43). 

= H Hi})  \r 

Equation  (44)  gives  the  volume  fraction. 

€ ij  - (T**v-  - CaAL) 

The  pressure  behind  the  flame  front  is  obtained  by  applying  a modified 
perfect  gas  relationship  to  the  total  gas  generated  from  every  cell. 

-?  « ( xtc  • » + p,  a,  %u) 

C Ml  - V»>*  - XTftJ) 


* •‘P-V^Z-r-  w"!*iS«^*T^^ 


where 

YT&  « SUn  of  gas  generated  in  each  cell  at  time  t2 

- distance  flame  front  has  progressed  from  t » y to  t^ 
v»u « total  volume  of  solids  behind  flame  front 
\ * propellant  covol  ime 

The  pressure  given  by  equation  (45)  applies  to  the  entire  region  behind 
the  flame  front.  This  pressure  is  averaged  with  the  pressure  from  the  preceding 
time  period  to  determine  the  flame  front  velocity.  The  procedure  is  repeated 
until  the  pressure  ir  equation  (45)  for  two  successive  passes  agrees  to  within  5X 
The  flame  front  calculation  procedure  has  been  programmed  and  the  details 
of  the  program  as  presented  in  Appendix  B. 
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Heat  Transfer  Models 

Since  the  analytical  model  considers  the  problem  in  three  distinct 
regions,  the  breech,  the  interior,  and  the  projectile  base  the  heat  transfer 
models  will  be  considered  in  the  same  manner.  The  approach  in  defining  the 
heat  rates  has  been  to  use  the  best  data  either  experimental  or  analytical. 

Breech  Region 

In  the  breech  region,  the  hot  gases  stagnate  and  as  such,  no  boundary 
layer  is  formed.  Although  there  may  be  recirculation  flows  in  this  region, 
analytical  description  would  be  extremely  difficult.  However,  there  is  data 
available  on  heating  in  the  projectile  case  region.  Reference  (12)  developed 
experimental  data  on  the  heating  rate  at  various  locations  in 
ammnition  case  region  or  a function  of  time.  Figure  5 presents  the  results 
from  reference  12.  To  obtain  the  heating  at  the  breech,  the  data  in  figure  5 
was  extrapolated  to  x * o and  then  replotted  in  figure  6.  The  peak  heating 
occurs  at  .25  cUsec  with  a value  of  3.7  Btu/ft^sec.  After  reaching  this 
peak  value  the  breech  heating  decreases  exponentially  with  time.  Rather  than 
trying  to  match  the  peak  heating,  the  approach  taken  was  to  satisfy  the  long 
term  heating  by  fitting  the  following  equation  to  the  data  in  figure  6. 

Using  i and  2 milsec  as  the  matching  times,  the  following  values  for  Ki  and 
K2  were  obtained: 

Ki  = 3.24  Htu/f t2sec 

K2  « “587.78  iec*1 

The  dotted  line  in  figure  6 gives  the  profile  of  the  fitted  curve.  To 
determine  the  effect,  of  breech  heating,  the  value  of  Kj  was  varied  trom  0 to  3.24 
on  a series  of  computer  runs.  The  results  showed  that  breech  heating  had  a 
negligible  effect  of  projectile  launch  velocity  and  peak  pressure. 


Interior  Region 

The  interior  region  refers  to  those  locations  not  considered  to  be  the 
breech  or  projectile  base.  In  the  interior  region,  the  boundary  layer  along 
the  wall  is  fairly  well  defined.  For  this  reason  the  approach  has  been  to 
use  the  standard  Reynolds  Analogy: 

» - £ 
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The  selection  of  the  friction  factor,  f,  is  somewhat  arbitrary  but 
values  between  .001  and  .006  have  been  used. 

Projectile  Base  Region 

Describing  the  heat  transfer  in  this  region  is  very  difficult  because 
the  normal  boundary  layer  approximations  are  not  satisfied.  In  this  region 
the  gas  velocity  at  the  wall  or  base  of  the  projectile  equals  the  projectile 
velocity  which  is  the  highest  velocity  in  the  problem.  In  addition,  the  level 
of  projectile  heating  (barrel  heat  transfer  at  projectile)  has  a profound 
effect  upon  the  peak  projectile  base  pressure  and  therefore,  the  launch  velocity. 

Initially,  the  approach  was  to  simply  use  the  Reynold-  Analogy  in  this 
region  with  the  gas  velocity  equal  to  the  projectile,  velocity.  However,  this 
approach  grossly  underestimates  the  peak  pressure  at  the  base.  Since  the 
pressure  increases  with  entropy  decrease  assuming  constant  gas  temperature, 
this  condition  is  valid  as  ljng  as  the  propellant  is  burning.  Entropy  can 
be  decreased  by  higher  heat  transfer  rates.  The  approach  taken  therefore, 
was  to  arbitrarily  increase  the  heat  transfer  coeff i<-ient . Factors  of  ten 
were  necessary  to  significantly  increase  the  base  pressure.  The  requirement 
of  a factor  of  this  size  seems  unreasonable.  In  addition  the  peak  heating 
rate  occurs  near  the  barrel  exit  which  does  not  agree  with  test  data. 


The  next  approach  was  to  use  the  heating  rate  based  upon  a shock 
tube  end  condition.  This  procedure  was  developed  in  reference  (13).  In 
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this  formulation  the  following  equation  was  used: 

U *■ 

In  this  representation  the  shape  of  the  heat  transfer  coefficient  as 
a function  of  time  is  correct.  However,  the  overall  levels  are  too  low. 
The  constant,  K,  can  vary  by  a factor  of  2 (reference  (13))  but  this  is 
still  Inadequate  to  predict  the  required  levels  of  heating. 

The  most  promising  procedure  is  the  one  described  in  reference  (3) . 
In  this  reference  the  flow  in  the  barrel  is  examined  qualitatively  and 
quantitatively.  The  reference  notes  that  away  from  the  projectile  the 
Reynolds  analogy  could  be  used.  However,  in  the  projectile  base  region 
the  momentum  thickness  goes  to  zero  because  the  gas  velocity  and  pro- 
jectile velocity  are  equal.  If  a boundary  layer  analysis  were  attempted 
this  would  yield  an  infinite  Nusselt  number. 

The  approach  taken  in  reference  (3)  was  to  use  the  usual  boundary 
layer  method  away  from  the  projectile  and  breech.  These  locations  are 
handled  using  approximate  solutions.  The  analysis  develops  the  integral 
momentum  equation  and  solves  this  equation  introducing  empirical  data 
when  needed.  The  projectile  heat  transfer  is  given  by: 

//„=  P* i 


imvnf  1 I'  • ■ * 
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where: 


A/*. 


Nusselt  Number  hcl/k 


* shear  law  constant  * .332 


*"  « Prandtl  Number 


Reyrolds  Number  based  upon  gas  velocity  and  L 


^ * Reynolds  Number  based  upon  momentum  thickness  at 
transition 


* Distance  from  projectile  to  breech 


The  transition  Reynolds  R.^  uken  „ 20„  Th#  hM{  tra„s£er  prediction 

method  of  reference  (3)  was  used  in  the  computer  model  because  it  provided  the 
right  profile  and  heat  transfer  rates  of  the  right  magnitude. 
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Numerical  Results 

Using  the  analysis  described  in  the  preceding  sections,  an  analytical 
model  computer  program  was  developed  and  number  of  cases  vere  examined. 

The  most  important  parameter  affecting  the  pressure  and  velocity  history  is 
the  surface  concentration  density  (SCD)  of  deterrent.  Figures  7 and  8 present 
the  pressure  and  velocity  histories  for  an  SCD  of  0.1  and  0.2  respectively. 

In  these  two  cases  it  is  assumed  that  the  flame  velocity  through  the  propellant 
bed  is  infinite.  The  surface  concentration  density  of  deterrent  refers  to  the 
distribution  of  deterrent  through  the  propellant  grain.  For  all  cases  using 
the  same  propellant,  the  total  amount  of  deterrent  is  the  same;  however,  for 
an  SCD  of  0.1  compared  to  0.2  the  deterrent  is  spread  to  a greater  grain 
depth.  This  means  that  the  effect  of  the  deterrent  is  less  over  a greater 
depth.  For  instance  the  SCD  - 0.1  deterred  flame  temperature  is  4615°R 
compared  to  3331  °R  for  an  SCD  *0.2.  From  figure  7 it  can  be  seen  that  the 
launch  velocity  for  an  SCD  * 0.1  is  2500  fps.  Launch  occurs  1.27  milsec  after 
ignition.  For  SCD  « 0.2,  (figure  8)  the  launch  velocity  is  1700  fps  at  1.69 
milsec.  For  both  these  runs  the  heat  transfer  rates  at  the  projectile  base 
were  defined  by  the  results  of  reference  (3).  Comparing  the  results  from 
figure  7 and  8 with  figure  1,  the  experimental  data,  it  can  be  seen  that 
considerable  differences  exist.  In  figure  7 the  launch  velocity  and  peak 
projectile  base  pressure  are  closer  to  the  experimental  results  than  figure  8. 
However,  the  launch  velocity  is  only  78%  of  the  measured  data  and  the  peak 
projectile  pressure  is  only  59%  of  the  corresponding  experimental  results. 

In  contrast  the  peak  breech  pressure  from  figure  7 is  183%  of  the  results  in 
figure  1.  It  is  felt  that  these  differences  can  be  attributed  to: 

(1)  projectile  heat  transfer 

(2)  flame  front  velocity 
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Based  upon  earlier  heat  transfer  studies  it  was  determined  that 
projectile  heat  transfer  had  a strong  influence  on  projectile  velocity 
and  projectile  base  pressure.  For  comparison  the  projectile  heat  rate  for 
an  SCD  of  0.1  and  0.2  are  presented  versus  travel  in  figure  9.  The  peak 
projectile  heating  is  14,500  Btu/ft^  sec  at  approximately  7 inches  of 
travel.  This  is  approximately  one  half  the  predicted  peak  from  reference  (3). 
in  addition,  the  peak  occurs  further  downstream  than  anticipated  by  refer- 
ence (3)  and  experimental  experience. 

The  large  variation  in  peak  breech  pressure  was  attributed  to  the 
h^gh  flame  speed.  To  examine  this  effect,  a finite  flame  program  was 
written  and  run.  The  analytical  basis  for  this  program  was  described  earlier. 
Table  1 presents  some  results  from  the  finite  flame  velocity  analysis. 

TABLE  L 


SCD  =0.1 

SCD  =0.2 

SCD  = o.; 

Max  Pressure 

414081 

14,720 

16,972 

boh ine  projectile-psia 

Flame  Travel  Time  - sec 

.595xl0-3 

.756xl0-3 

.731x10 

SCU=0. 1 but  1/2  burning  rate  to  simulate  ignition 
delay 

Figure  10  presents  the  volume  fraction,  r. , distribution  versus  distance 
for  the  same  cases  on  Table  1. 

The  pressure  and  velocity  history  using  a finite  flame  velocity  is 
given  in  figure  ii  for  SCD=0.2.  The  launch  velocity  is  increased  by  400  fps 
and  the  projectile  base  pressure  is  increased  slightly.  The  breech  pressure 
is  reduced  significantly  because  the  finite  flame  velocity  analysis  sets  Ihu 
volume  fraction  at  the  breech  equal  to  zero.  Figure  12  presents  the  same 
results  for  an  SCD=Q.l.  The  launch  velocity  in  figure  12  is  3800  fps. 

The  peak  projectile  base  pressure  is  64,500  psia.  This  represents  a sub- 
stantial change  over  the  results  in  figure  7.  The  flat  response  of  the 
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breech  pressure  curve  is  again  attributed  to  the  zero  volume  fraction  at 
the  breech. 

To  eliminate  the  fiat  response  at  the  breech  when  using  the  finite 
flame  criteria,  it  was  decided  to  set  the  volume  fraction  at  breech 
equal  to  the  first  station  in  the  flame  front  analysis.  This  change 
Increases  the  breech  pressure  significantly.  For  an  SCD  of  0.1,  the 
results  are  not  reasonable.  To  moderate  this  effect,  the  one  half  burn- 
ing rate  was  used  in  the  flame  front  analysis  (see  Table  1).  With  the 
breech  volume  fraction  equal  to  zero,  the  peak  breech  pressure  was 
approximately  17*000  psia.  The  peak  projectile  bore  pressure  and  launch 
velocity  were  57,500  psia  and  3640  fps  respectively.  Figure  13  presents 
the  complete  pressure  and  velocity  history  for  the  non  zero  breech  volume 
fraction.  In  this  figure  it  can  be  seen  that  the  peak  breech  pressure  is 
130,000  psia.  This  is  approximately  twice  the  experimental  results.  For 
this  case  the  starting  breech  volume  fraction  was  0.574.  A lower  value 
would  reduce  this  pressure. 

Regarding  the  choice  of  heat  transfer  model,  a number  of  runs  were 
made  using  the  usual  Reynolds  analogy.  The  results  from  these  runs  yielded 
launch  velocities  and  pressure  at  the  projectile  well  below  test  data. 

Some  parametric  runs  were  made  to  determine  the  necessary  increase  in  heat 
transfer  required  to  yield  agreement  with  test.  These  results  indicated 
a factor  of  10  or  20  which  does  not  seem  reasonable.  In  addition  the 
general  shape  of  the  profile  predicted  by  the  Reynolds  analogy  is  in- 
correct. The  Reynold's  heating  peaks  with  peak  velocity.  This  means  that 
peak  heating  would  occur  near  launch.  For  this  reason,  the  decision  was 
made  to  use  the  heat  transfer  model  described  in  inference  (3)  exclusively. 
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Conclusions  and  Recommendations 

Computer  models  have  been  used  to  describe  the  flow  field  in  small  ar 
and  also  to  describe  the  propagation  of  l lame  front  through  a porous 
propellant  bed.  The  model  determines  the  pressure,  velocity,  propellant 
volume  fraction,  and  temperature  versus  time  and  position. 

Numerical  results  have  shown  the  dependence  of  projectile  velocity 
and  barrel  pressure  history  upon  deterrent  concentration  density  and 
flame  velocity.  It  has  also  been  determined  that  the  heat  transfer  to  the 
base  of  the  projectile  also  has  a significant  effect  upon  the  launch 
condition. 

Unfortunately,  it  was  not  possible  to  simulate  the  experimental  data 
presented  in  figure  1.  However,  with  some  additional  effort,  the  author 
feels  that  the  results  of  figure  1 can  be  duplicated.  Toward  this  goal, 
tie  following  recommendations  are  made: 

(1)  Vary  breech  volume  fraction  to  simulate  breech  pressure. 

(2)  Parametric  study  of  projectile  base  heat  transfer  to  obtain  a steeper 
pressure  profile. 

( 1)  Parametric  study  of  projectile  resistive  force. 

(4)  Development  of  analytical  equations  for  flow  at  the  base  of  the 
projectile. 
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COMPUTER  PROGRAM 
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The  computer  program  implementing  the  numerical  procedures  describer  in 
preceding  sections  composed  of  a main  program  and  10  subprograms.  The  main 
program  controls  the  flow  of  the  calculation.  The  subprogram  performs  certain 
specific  functions  during  the  calculation.  The  name  and  function  of  each  sub- 
program is  given  below; 


SUBPROGRAM 

DIME 


START 


INTER 

BREACH 

BASEP 

OUTPUT 

DIAGO 


RFORCE 

PROP 


FUNCTION 

calculates  the  surface  area,  Sp,  and  volume, 

Vp,  of  propellant  grain  as  a function  of  dis- 
tance normal  to  the  surface. 

generates  the  pressure,  entropy  time  history 
along  the  first  Q wave  from  the  starting  point 
to  the  breech.  The  required  number  of  points 
is  specified  in  the  input  and  program  divides 
the  distance  from  projectile  base  to  breech  into 
equal  increments. 

calculates  the  interior  points  using  value  at 
(3)  and  (4)  (figure  2)  as  known 

calculates  the  breech  points 

calculates  the  projectile  base  points 

points  the  required  output  information  in  the 
specified  format 

a diagnostic  program  which  prints  certain  informa- 
tion if  the  calculations  in  INTER,  BkEACH,  and 
BASEP  do  not  converge  within  10  tries 

calculates  the  engraving  force  as  to  function 
of  X based  upon  tabular  input  data 

calculates  the  specific  heat  ratio,  gas  constant 
and  specific  heat  at  constant  volume  as  a function 
of  the  amount  of  deterred  and  ncn-detcrrcd  pro- 
pellant burned.  The  program  is  called  only  after 
the  deterred  layer  is  burned  through 
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TAREA  This  is  a subjnnction  which 

calculates  the:  tube  cross 
sectional  area  as  a function 
of  X based  on  tabulated  inputs. 


The  main  program  maintains  control  of  the  overall  program  and  stores 
the  calculated  values  of  various  parameters  (i.e.,  gas  velocity,  propellant 
velocity,  pressure,  etc.)  This  data  is  3tored  in  a double  subscripted  (i,j) 
array.  When  i equals  j,  the  points  correspond  to  the  base  of  the  projectile. 
Along  a given  Q wave  the  j subscript  is  constant.  For  a fixed  value  of  j, 
the  number  of  i values  which  determines  the  grid  size  is  an  input  parameter. 
Figure  5 present'*  in  abbreviated  flow  diagram  for  the  main  program. 
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X(I.J) 

TIME (I, J) 

VEL(I,J) 

SS(I,J) 

I’RESS  ( I,  J) 

ICMP(I.J) 

DBWS(I,J) 

VELP(1,J) 

ENTR(1,J) 

V0LFRA(1,J) 

AJ.(I , J) 


X,  location  of  grid  points 

t , time  of  grid  points 

gas  velocity  at  grid  points 

speed  of  sound  at  grid  points 

pressure  at  grid  points 

gas  temperature  at  grid  points 

gas  density  at  grid  points 

propellant  velocity  at  grid  points 

entropy  of  gas  at  grid  points 

volume  fraction  at  grid  points 

distance  burned  normal  to  propellant 
gram  surface  at  grid  points 


The  following  are  used  in  the  main  program  and  subroutine  INTER  (see 
Figure  2) 


X3,  X.4,  Xb 
T3,  T4,  T6 
UG3,  UG4,  UG6 
A3,  A4,  A6 
PSI3,  PS] 4,  PS16 
AREA3,  AREA4,  AREAS 
TG3,  TG4,  TG6 
RH03,  RH04,  RHOG 


Spatial  location  of  points  3,  4 and  6 

time  location  of  points  3,  4 and  6 

gas  velocity  at.  points  3,  4,  and  6 

speed  of  sound  at  points  3,  4 and  6 

gas  pressure  at  points  3,  4 and  (- 

tube  cross  sectional  area  at  points  3.  4 and  6 

gas  temperature  at  points  3,  4 and  6 

gas  density  at  points  3,  4 and  6 


UPS,  UP4,  UP6 
S3,  S4,  S6 
EPL3 , EPU,  EPL6 

L3,  L4,  LS,  1,6 


propellant  velocity  at  points  3,  4 and  6 

entropy  at  points  3,  4 and  6 

propellant  volume  fraction  at  points 
3,  4 and  6 

distance  burned  normal  to  propellant 
grain  at  points  3,  4,  S and  6 


The  following  are  used  in  the  main  program  and  subroutines  BREACH 


and  START  (see  figure  4) 
X12 

Til,  T12,  T13 
UG11,  UG12,  UG13 
All,  A12,  A13 
PSI1,  PSI2,  PS I 3 
TGI  1 , TGI 2,  TGI 3 
PJ10U,  RH012,  WJ013 
SJ1,  S12,  S13 
EPL11,  EPL12,  KPL13 

AKEA11,  AREA1 2,  AREA 13 


spatial  location  of  point  2 

time  location  of  point  1,  2 and  3 

gas  velocity  at  points  1,  2 and  3 

speed  of  sound  at  points  1,  2 and  3 

gas  pressure  at  points  1,  2 and  3 

gas  temperature  at  points  1,  2 and  3 

gas  density  at  points  1,  2 and  3 

entropy  at  points  1,  2 and  3 

propellant  volume  fraction  at  points 
1,  2 and  3 

tube  cross  sectional  area  at  points 
1 , 2 and  3 


in,  ns 


distance  burned  normal  to  propellant 
grain  at  points  1 and  3 


UP12 


propellant  velocity  at  point' 2 


The  following  are  used  in  the  main  program  and  subroutine  P-ASEP  (see  figure  3) 


X21 , 

X22, 

X24 

T-l  % 

'T*"> 

1 *.*• 

spatial  location  of  points  1,  2 and  4 
time  location  of  points  1,  2 and  4 
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UG21,  UG22,  UG24 
A?l,  A22,  A24 
PS121,  PSI22,  PS I 24 
TG21,  TG22,  TG24 
HI. 021,  RH022,  RH024 
S.?1 , S22,  S24 
FPL21,  CPL22,  HPL24 


AREA21,  AREA22,  AREA24 


L21,  L24 


NSTART 


gas  velocity  at  points  1,  2 and  4 

speed  of  sound  at  points  1,  2 and  4 

gas  press:. re  at  points  1,  2 and  4 

gas  temperature  at  point.  3,  2 and  4 

gas  density  at  points  1 , 2 and  4 

entropy  at  points  1,  2 and  4 

propellant  volume  fraction  at  points 
1,  2 and  4 

tube  cross  sectional  area  at 
points  1,  2 and  4 

distance  burned  normal  to  the 
propellant  surface  at  points 
1 and  4 

the  value  of  the  I index  for  the 
Q wave 
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subroutines: 


Name : 


J nput : 


Output : 


Method: 


Name : 


Input : 


INTER 


Purpose : To  calculate  the  characteristics  and  thermodynamic  properties 


at  any  interior  point  in  the  mesh 


X3,  X4,  T3,  T4,  UG3,  UG4,  A3,  A4,  PSI3,  PSI4,  TG3,  TGi,  RH03 


Rl!04;  UP3,  UP4,  S3,  S4,  EPL.3,  EPL4,  AREAS , ARF.A4 , L3,  L4 


X6,  T6,  UG6,  A6,  PSI&,  TOO,  TH06,  UPC,  S6,  EPI.6,  AREA 6, 


An  iterative  shcerae  is  used  until  convergence  to  within 


a given  tolerance  is  obtained.  Normally,  the  convergence. 


check  is  made  on  the  calculated  speed  of  sound  and  gas  velocity. 


However,  for  gas  velocities  below  300  ft/sec  only  the  speed 


of  sound  is  checked.  This  procedure  was  necessary  because 


below  300  ft/sec  the  accuracy  requirements  on  p and  Q for 


calculating  gas  velocity  become  excessive.  At  300  ft/sec, 


a one  percent  error  in  the  speed  of  sound  can  be  translated 


into  approximate  a 10  percent  in  the  gas  velocity.  At  higher 


gas  velocities,  the  percentage  tolerance  on  gas  velocity  and 


speed  of  sound  are  the  same. 


BREACH 


Purpose:  To  calculate  the  characteristics  and  thermodynamic  propei tie.': 


at  the  loc.u  ion  X=0 


XI 2,  Til,  112,  UG11 , UG12,  All,  A12,  PSU1,  PS112,  TGi  1 , TGI 2. 


RiiOll,  RH012 , Sll,  S12,  LP1.11,  r.PL12,  AREA  11 , AREA! 2,  Ul, 


II I AM  MM  I 


Output: 

Method: 


T13,  UG13 , A13,  PSI13,  TG13,  RH013,  S13,  EPL13,  1.15 
An  iterative  scheme  is  used  until  convergence  to  within  a 
given  tolerance  is  obtained. 


Name: 


Purpose: 


Input : 


BASliP 

To  calculate  the  characteristics  .and  thermodynamic  propevties 
at  the  base  of  the  projectile  and  to  locate  the  base  of  the 
projectile 

X21,  X22,  T21 , T22,  UG21,  UG22,  A21,  A22,  PSI21,  PSI24, 

TG22,  RH021 , RH022,  S21,  S22,  EPL21,  EPL22,  AREA21,  AREA22, 
L21 , NSTART 


*'*'"  * * * "***’ 
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Output : X24,  T24,  UG24,  PS124,  TG24,  UH024,  S24,  EPL24,  ARE A24 , 1,24 

Method:  An  iterative  scheme  As  used  until  convergence  to  within  a given 

tolerance  is  obtained 


Name : 


Purpose: 


Input : 


Out  put : 


Method : 


DIME  (A,  R,  SURA,  SURB,  VOLA,  VOUi) 

To  calculate  the  average  propellant  grain  surface  area  and 
volume 

A,  B corresponding  to  the  distance  burned  normal  to  the  surface 
( jL  ) at  two  points 

SURA,  SURE  surface  of  gra^n  at  respective  points 
VOLA,  VOLB  volume  of  grain  at  respective  points 
The  following  equations  are  used 
Volume  “ '2.  rr  C ^ ~ 

+ ( ( g-S?;}  +±( 

3>  XT 

Surface  area  - •ztt  ( 1«  C Rc') 

»■  ..  « * * 


*-  (6  - F-yY) 


name : 


OUTPUT 


Purpose:  To  print  the  calculated  results 

Input:  Calculated  velocities  and  theimodynamic  properties  at  all 

of  the  grid  points 

Output : Printed  results  including  values  of  P and  Q at  each  grid  point 

Method:  Format  of  printed  results  specified  in  subrout ’no 


Name : 


DIAGO  (LA,  BUG,  AA) 


Purpose:  A diagonostic  tool  to  print  results  calculated  prior  to 


failure  to  converge  event 
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input : 


Output 


Method : 


Name: 

Purpose: 


Input : 

Output: 

Method: 


Name: 

Purpose: 


Input: 


Output : 
Method: 


LA  indication  of  which  subroutine  failed  LA-1  INTER 

LA-2  BASHP 
LA-3  BREACH 


BUG  last  calculated  gas  velocity 
AA  last  calculated  speed  of  sound 

Printed  output  up  to  and  including  current  value  of  thermo- 
dynamic variables 

'ormat  used  is  specified  in  the  subroutine 


RF0RCB(A,B) 

Calculates  the  engraving  force  or  a function  of  distance 
from  initial  projectile  location 

A,  location  of  the  base  of  the  projectile 

B,  resistance  force  in  lbf 

Resistance  force  tabulated  or  function  of  distance  based  upon 
data  from  Report  R-2066,  ’’Study  of  the  Pressure  Distribution 
Behind  the  M193  Projectile  Vhen  Fired  in  the  M16  Rifle  Barrel" 


PROP (A) 

Calculates  the  CV,  GAMMA,  and  RGAS  as  a function  of  the 
percentage  deterred  and  non-deterred  propellant  burned. 
Subroutine  is  used  only  when  distance  burned  normal  to  the 
surface  is  greater  than  the  deterred  thickness  of  C1  ■•00084  ft 
A distance  burned  normal  to  the  propellant  surface  at  a given 
mesh  point 
CV,  GAMMA,  RGAS 

Calculation  assumes  the  gas  generated  to  be  ideal  and  that 
properties  are  based  upon  mass  fraction  of  components  present 


40 


Name : 


START 


Purpose:  Calculates  the  thermodynamic  properties  along  the  first  Q wave 

at  specified  intervals 

Input:  Same  as  BRRACll  but  here  they  refer  to  data  along  first  Q wave 

Output : Same  as  BRHACH  but  here  they  refer  to  data  along  first  Q wave 

Method:  The  subroutine  is  called  a specified  number  of  time  (i.o.  an 

input  parameter).  Bach  time  it  calculates  the  therraodyne^ic 
properties  along  the  first  Q wave.  The  results  are  equally 
spaced  between  the  projectile  initial  location  and  the  breech. 
The  gas  velocity  and  heat  transfer  along  this  wave  are  zero 
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Name:  TARGA(ALOC) 

Purpose:  Calculates  the  tube  cross  sectional  area  as  a function  of 

distance  from  the  breech 
Input:  ALOC  location  of  the  grid  point 

2 

Output:  TAREA  tube  cross  sectional  in  ft 

Method : Interpolates  data  from  table  which  is  supplied 


I * ’ 
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C THIS  PROGRAM  CALCULATES  THE  PRESSURE . TEMPERATURE  AND  VELOCITY 
C TIME  HISTORIES  BASED  UPON  A GIVEN  PRCPELLANT  CHARACTERISTIC 
C THIS  IS  THE  MAIN  PROGRAM 

REAL  L5fL6fl3,L4,L2ifL24,LllfL13 

COMMGN/C/XI 99*99) (TIME 199,99) ,VE LI 99,99 ), SSI  99, 99), PRESS 1 99*99 ) ,TE 
$MPI99,99),DEN$( 99*99) , VE IP (99,99 ) ,ENTR<99, 99 ) , VOLFRA 1 99, 99) ,AL(99, 
799) 

COMMON/CCM1/X3, X4,X6,T3,T4,T6,UG3,UG4,UG6, A3,A4,A6,PSI3,PSI4,PSI6, 
STG3»TG%»TG6,RH03»RH04,RHC6»UP3,UP4»UP6,S3,S4,$6«EPL3»EPL4*EPL6*ARE 
>A3, AREA4, AREA6,L3«L4«L5, L6 

COMMCN/CCM2/X12,TLl,T12, T13,UG11 ,UG12,UG13t All, A12,A13*PS11 l,PSI 12 
$,PSI 13,TG11 «TG12*TG13,RHC11*RH012,RH013,S11, S12,S13,EPL11,EPL12,EP 
*L13,AREA11 , AREA12, AREA13,Lll,L13,UP12 
COMMGN/CGM3/X21 «X22,X24, T21,T22, T24,UG21 »UG22,UG24, A21,A22« A24,P$1 
121 »PSI22»PSI24*TG2l»TG22tTG24tRHC2l«RH022*RH024»S21f S22*S24»EPL21t 
)EPL22»EPL24»AREA21* AREA22, AREA24»L21»L24»NSTART 

COMMCN/SAME/RGAS,PSI 1 «GAMMA,DZERC,PIE,6EE, PODPOT,CV, BASE,GZERO, WT • 
>TOL* AN 

COPMCN/STAR/KSTART 
COMMON/OUT/MAXJtMAXI 
COMMCN/EXIT /XLIM»TMUZE 
COMMON /6HEAT/F ACT, TWALL,FRICT 

COMMON/PARM/PAR1N1 • PAR IN2 1 PAR IN3 , PARBA1 , PARBA2 . PARBA3, PARBR 1 , PARBR 
<2t PARBR3 

COMMON/PCWOER/CV1 tCV2*  RGASl* RGAS2»GAMlt GAM2« OLAY 

NAMEL I ST/PAR/PAR  INI, PAR I N2, PAR IN 3,PARBA1,PAR8A2,PARBA3, PARBR 1»PARB 
SR2.PARBR3 

NAME LI  ST/ IN POT/ BASE  ,WT,TCL,XLIM, PSI 1, DZ£RO«TZERQ, XZERO,FRACl, ALLl, 
IFACTtTNALL* AN«FRICT , I , J 

NAMELI ST/PRCPEL/PP1 »PP2*CV1,CV2»  BEEl,B£E2»RGASl»RGAS2,DLAY»GAMi,GA 
•M2 

999  REA015, INPOT) 

REAOlSf PAR) 

REA0I5, PROPEL I 
WRITEI6, PROPEL) 

KRITEI6, INPCT ) 

MRITEI6«PAR ) 

00  1 II»l,99 
DO  l JJ*1,99 
SS  1 1 1 , J J ) -0  .0 
1 CONTINUE 

TMLZE-1000. 

P00PCT«PP1 

CV»CVl 

GAMMA-GAMl 

BEE-BEE1 

RGAS-AGASl 

PIE-3.141593 

GZERC-32.2 

TATM-530. 

IF ( I .GT.O )GC  TO  20 
C CALCULATION  CF  FIRST  TWO  POINTS 
XI 1, 1 ) -XZERC 
X<2,1)«0.0 
TIME! ltl )*0.0 

SSI 1 1 1 )»SQRTIGAMMA*RGAS*TZERO) 

VOLFRAI 1 , l ) «FRAC1 

VEL( If  1 )*0.C 

VELP(l»i)*0.0 

PRESS ( It  1 )«P$I 1 

TEMPI!,  U-TZERC 

OENSI 1, 1 )*PSI 1/RGAS/TZERC 

ENTRI 1 , 1 )*( RGAS/778.+CV I ♦ALOGITZERO/TATMJ-RGAS/778.PALOGIPSI 1/2075 


Mwm  '->r*A.»-».it«. 


#•) 

ENTRU,i)*EKTR<l,l;  *77€WRGAS/GAPPA 
AL(ltl)*ALLl 

IF(ALLl.GT.CLAY)  BEE*BEE 1 
IF( ALL1.GT.CLAY)  P0DP0T-FP2 
IFCALLi.GT.CLAYKAU  PRCPIAUl) 

KSTART-J 

DO  3 1*1, KSTART 

K*m 

X12*X< 1*1) 

T12»TIHE(I*1) 

EPtl2*V0LFRA< 1,1) 

A12*SS<Itl> 

UG12*VELU,1) 

UP12*VELP<I,1) 

psii2»PRESS(i*n 

TG12«TERPU»1) 

RH012*DEN$< 1,1) 

S12*ENTR( I , 1) 

Lll*AL( 1,1) 

CALL  START 
PRPSI-PSI13/144. 

WR!TE(6,5002) A13*L13»PRPS1,TG13,RF013«S13,EPL13,T13 
X( K, l )«X12 
TIPE(K,1)»T13 
VOLFRA(K, 1 ) “EPL13 
SS(K, 1 )*A13 
VEL(K,1)*UG13 
VELP( K,1 )*UG13 
PRESS (K« 1 )*PSI 1 3 
TEMP(K,i)»TG13 
DENS ( K , 1 )*RH013 
ENTR ( K, 1 )*S13 
AL ( K , 1 )*L13 
3 CONTINUE 
NCHECK-0 
JLAST«1 

5 DO  50  1*2,80 

IF(NCHECK*GT.O)JLAST«JLASTU 

NCHECK-0 

JQ* JLAST+1 

00  A9  JJ*J0,I 

J-JJ 

IF(JJ.EG.I)  GO  TO  55 
X4*X ( I-i » J ) 

X3*X(I, J-l) 

T4*TINE( I-i »J) 

T3*TINE( I, J-l) 

UG«*VEL(I-1,J) 

UG3*VELU*J-1I 
A4«SS(I-lt J) 

A3*SS ( I ♦ J-l ) 

PSI**PRESS<  I-i,J) 

PSI3«PRESSC I, J-l) 

TG4*TENP( 1-1, J) 

TG3-TENPU,  J-l) 

RHC4*0ENS(I-l,J) 

RHC3*0ENS( I , J-l ) 


UP*»V€LPU-l,  J) 
UP3*VELP< I , J-i) 
S*»ENTRU-1,J> 
S3*£NTR { I , J-l ) 
£PLA*VCLFRA(I-i , J) 


1 “‘  ,,u 


I 

V 

f 


I 


I 


1 


i’ 

b 

iK 


f 


m 


o 

0 

0 

0 

0 

0 

D 


EPL3*VCLFRA!I,J-1) 

L4*AL i 1-1 » J ) 

L3*AL ! I t J-l ) 

T00L*!L3+L4)/2« 

IF (TCCL«GT«  CLAY ) GO  TO  1G 
GO  TO  14 
10  BEE*BEE2 
P0DP0T-PP2 
CALL  PROP(TCOL) 

14  AREA3-TAREAIX3) 

AREA4*TAREA!X4) 

WRITE!6t5003)X3,X4,T3,T4,UG3.UG4.A3.A4,PSI3,PSI4,TG3.T64#RH03,RM04 
$»UP3tUP4,£PL3fEPL4*L3*L4,AREA3fAREA4 
CALL  INTER 
RGAS*RGAS1 
BEE-BEE1 
GAMMA-GAM1 
CV«CV1 
P00PCT*PP1 
PRPSI-PSI6/144. 

WRITE(6,5002)A6,L6,PRPSI,TC6,RH06,S6tEPL6,T6tX6,UG6,UP6 
5002  FORMAT !1X#11E11«3) 

X(I, J)»X6 


I 


n 

L! 


; i 

i 

t J 


TIMEU# J»*T6 
VEL ( I « J )*UG6 
PRESS! I » J)*PSI6 
TEMP! I » J)*TG6 
OENS! I » J )*RH06 
VELP! I * J)*UP6 
ENTR! I » J)*S6 
VOLFRA! I f J)*EPL6 
SS< I * J )*A6 
ALU*  J)*L6 

IF(X3.LE.0.C)  GO  TO  40 
IFCT6.GE.TMLZE)  GO  TO  48 
GO  TC  49 
55  L* JJ-1 
NSTART*I~i 
K«l 

X2l*X(NSTART»L) 


X22*X(K*L) 


T2l«TIME(NSTART,L) 
T22»TIME(K,L) 
UG2I*VEL(NSTART,L ) 
UG22*VEL (K»  L> 
A2l«SS ( NSTART »L ) 


A22«SS(K*L) 
P$I21*PRE5$ (NSTART* L) 
PSI22*PRESS!K,L) 
TG21»TEMP!NSTART,L> 
TG22*TEMP(K  »L ) 


RHC2i-CEN5t  NSTART* L } 


RHC22»DENS< K.L) 

S21»ENTR!NSTART»L) 

$22«ENTR!K,L) 

EPL2l*VCLFRA<NSTARTtL) 

EPL22-VCLFRA(K,L) 

L2l*AL (NSTART*L ) 

IF ( L21 «GT.DlAY I BEE-BEE2 
IF (L21«GT«0LAY)  P00PCT*PP2 
IF(L21.6T«0LAY)  CALL  PR0P!L21) 
AREA2l*TAREA<  X21 ) 

AREA22«TAREA! X22) 


CALL  8ASEP 
NSTART*I 
P0DP0T-PP1 
CV=CV1 
GAMMA*GAMi 
8EE-BES1 
RGAS=RGAS1 
PR PSI*PSI 24/144* 

WRITE <6, 5001 )A24,L24,PRPSI,TG24,RH024,$24, EPL24.T24, X24.UG24 
5001  FORMAT (1X*10E11*3) 

X<NSTART,J)*X24 
TIME<NSTART,J)*T24 
VEL:NSTART*J»*UG24 
V6LP(NSTART,J)*UG24 
$S(NSTART*J)*A24 
PRESS(NSTART*J)*P5I24 
TEMP (NST ART  * J)*TG24 
DENS ( NStART  * J )*RH024 
ENTR(NSTART,J)»S24 
VOLFRAI NSTART * J )*EPL24 
AL(NSTART,J)*L24 
GO  TO  49 
40  NCHECK*1 
II*I 

X12«X( lit J) 

Tli*T13 

T12=T I ME ( 1 1 * J ) 

UG1 1*0*0 
UG12*VEL ( t I » J) 

X*1=0.0 

A11-A13 

4i2*S$lII,J) 

PSI11*PSI 13 
PSU2*PRESS<II,  J) 

TG11*TG13 
TG12*TEMP ( 1 1 * J ) 

RHC11»RHC13 
RH012=0ENS( Ilf J) 

SI  1=  SI  3 

S12=ENTR (II  ,J) 

EPLU*EP'.13 
EP112»VCLFRA<II ,J> 

UP12*VELP( i I*J) 

LI  i = H 3 

IFU11.GUDLAY)  8E6*8EE2 
IFail.GT.OLAY)  P0DP0T*PP2 
I F ( L 1 1 *GT* OLAY ) CALL  PROPlLil) 

AREA11*TAR£A(X11) 

AREA12=TAREA<X12) 

CALL  BREACH 

RGAS*RGAS1 

8EE*8£E1 

r.Aku»*r.«u) 

CV*CV1 
P0GPCT*PP1 
PRPSI*PSI 13/144* 

WP. I TEI  6*5002 ) A1 3f  L13*  PRPSI  * TG13*  RF013* S13*  EPL 13*T13*  X13 
L = I ♦! 

K*  J 

AL(LfK)*L13 
TI  ME  (LfK)s»T  13 
SS1L,K)*A13 
PRESS! L»K)*PSI13 


»-  - TIJ*'  "V  I 


' r'*'*V-V/5P  ^ R^***1'  l 


.,  v»  ./'*'.'»*  **** 


J J 


1 i 


49 
48 

50 
56 


TEMP(LtK)*TG13 
0ENSU,K)*RH013 
VEUL,K>*0.C 
VELP(L,K)*C.O 
ENTR<L,K)*S13 
X<L,K)*0.0 
VGLFRA'L,K) *EPL13 
MAXJ*JJ 

IF(T13.G6.TMUZE)  GO  TO  46 

CONTINUE 

MAXIM 

IFtT13.GE.TMUZE>  GO  TO  56 
CONTINUE 
CALL  OUTPUT 
MAXI*MAXI«-1 
00  51  I*MAX J,MAXI 


00  51  I*MAX J,MAXI 

WRITE (7,5003) X(I ,MAXJ) ,TIME( l, MAXJ) , VOLFRA ( I , MAX J ), SS ( It 
I»MAXJ)»VELP(I» MAXJ ) t PRESS ( I »HAX J ),DENS( I , MAXJ ) , ENTR ( l ,M 
MAXJ) « TEMPI I, MAXJ) 

CODMITIICIA  Kl 


51 

20 


»TlMEU,l).VOLFRAU,l),SS(I,i>,VEUI.l>,V 
1*1) ,ENTR< I ,1 ) , AL( I, 1), TEMPI I, 1) 

) ,TIME (I »1 )* VOLFRA ( 1*1) *SS( 1 , 1 ) , VEL ( 1 , 1 ) , 
(1*1  )*  ENTR( I* 1 ) *AL ( 1*1 ) *TEMP( 1*1) 


SI 

#MAX  J 

5003  FORMAT ( 3E20*5) 

CONTINUE 
GO  TO  999 
NSTART*1 
KSTART*J 
JLA$T*1 
NCHECK-0 
KK* J4 1 

00  52  1*1, KK 
R£A0(5,5003)XU,1),TIME 
t)  , PRESS! 1*1  ), DENS (1,1) • 

WRITE<6,5003)X(I,1)-TIM 

#1 ) * PRESS! 1,1) , OENSI 
52  CONTINUE 

T13*TIME(KK  »1 ) 

A13*SSIKK, 1 ) 

PSI13*PRESS (KKt 1 ) 

TG13*TEMP{ KK* 1 ) 

RHC13*DENS(KK,1) 

S13*ENTR(KK,1) 

EPL13*VCLFRA(KK,1) 

L13*AL(KK,i) 

GO  TC  5 
ENO 

SUBROUTINE  8ASEP 
SUBROUTINE  BASER 

REAL  L5.L6  . L3,L4.L2l,L24,LU,L13 
THIS  SUBROUTINE  CALCULATES  PCINTS  AT  THE  BASE  OF  THE 
COMMON /COM 3 /X 21 ? X22»X24» T2 i * T22, T 24, UG21 »UG22* UG 
S21,PSI22,PSI 
JEPL22, EPL24 . 


. MAXJ ) , V£L( 
MAXJ), AUI, 


VELPI 1,1 
VELPII, 


CALCULATES  PCINTS  AT  THE  BASE  OF  THE  PRO 
/X21 ?X22»X24| T2I»T22*T 24,UG21 »UG22»UG24, A 
124,TG2l,TG22  ,TC-24, RHC21  ,RH022,RH024,S21, 

AREA21 , AREA22, AREA24 , L21,L24, N START 

/RGAS.PSIl. GAMMA. 02FRC. PIE. BEE, P00P0T,CV, BASE, GZERO,WT, 


PROJECTILE 

A21,A22,A24,PS! 
S22,S24, EPL21, 


) EPL22 , EPL24  AREA21 , AREA22, AREA24 , L21,L24, N START 
COMMON/SAME /RGAS, PS I 1, GAMMA, OZERC ,PIE, 

>TOL,AN 

pnuumucyrr  /vt  tM  Tkiuc 
y*ur  run/UAi  I / A W l T f I PULL 

COMMON/ CHEAT/ FACT »TWALL,FRICT 

COMMCN/PARM/PAR INI , PAR  IN2s  PAR IN3 , PARBA1 , PARBA2 , PARBA3, PARBR 1 »PARBR 
$2 , PARBR3 

OATA  RHCP,REC,0ELTT/3, 1085,0*9,0 ,000005/ 

Al*SCiRT  (GAMMA*RGAS*530, ) 

XTRAV*(X21-0.125)*12. 

WRITE <6, 5000 )POOPOT,CV,BEE,L21,XTRAV 
ITRY-1 

NSTART-NSTARTU 
5 PSI24-PS121 


■ ■>!■»  JIJUfH  WW7  1 


EPL24*EPL21 
A24*A21 

P22*2.*A22/ (GAMMA-1. )*LG22 
RHC24*RHC2i 
UG24-UG21 
TG24*TG21 

ABASE* BASE ♦ (PSI 24+PSI2 1 ) 

10  SL0PEP*<UG22*A22HlG244A2A)/2. 

CALL  RFCRCE ( X21 ,R1 ) 

C*X?1-X224$LCP£P*T22 
WORK* ( A8A$E/2*-Rl )*GZ£RC*T21/WT 
C»C-W0RK*T21/2.-T21*IUC21-WCRK) 

B*UG21-WCRK-SL0PEP 
A*(A8ASE/2.-Rl)*GZER0/WT/2. 

5000  FORMAT ( IX, • ♦*****♦«  t //.IX, 5E20. 5 ) 

T24* (-B-SQRTI B*8-4**A*C ) )/2./A 
X24*X22tSL0PEP*(T24-T22) 

IHX24.GE.XLIM)  GO  TO  75 

11  AREA24*TAREA( X24) 

AREA*! AREA24+AREA21 )/2. 

UGMH* < UG21+LG24 ) /2. 

PSIMN*(PSI244PSI21)/2. 

BR*BEE*( (PSIMN/144. )**AN  )/12. 

0P0X*0*0 
DL*BR*(T24-T21 ) 

L24*L21+DL 

IF (L24.GT.0. 000625 )L24*0.0C0625 
CALL  0IME(L21,L24,SURF21  , SLRF24, VCL21, V0L24 ) 
IF(VOL21.LE.O.OO)SURF21«C.OO 
IF(V0L2A.LE.0.00)$URF24»C.QC 
0BAR*0ZER0-L24-L2l 
IF (DBAR. LE. 0.000652 )0BAR*0.00 
0BAR*0BAR4*3 
F0R24»-PIE*CBAR*0P0X/6. 

EPLMN* (EPL21+EPL24 ) / 2. 

VOLMN* (V0L2 1+V0L24) /2. 

SURFMN* ( SURF21+SURF24) /2 . 

IF(VOLMN.LE *0*00 ) GO  TC  12 
ETAMN*EPLMN/VOLMN 
GO  TG  13 

12  ETAMN*0.00 

13  CONTINUE 

AMGMN«ETAMN*SURFMN*BR*RHCP*AREA 
DUPDX* ( UG24-UG2 1 ) / ( X24-X  21 ) 
0EPLCT*EPL21Ml.-EPL21>*<RF024-RK2i)/RHG2l-(EPl2im.-EPL21)*RHO2 
#1/RH0P )*AMGMN*<  T24-T21 )/RHC2l/AREA 
0EPL0T*DEPLCT/(T24-T21) 

14  CONTINUE 

FD24*EPLMN*F0R24/V0LMN 
RHOMN* < RHG244RHC21 ) /2. 

TGWN* ( TG24MG21 ) / 2. 

01  A* SORT  (4. *AREA/PIE ) 

AMN* t A22  + A24 ) /2 • 

TAW=TGMN*(1.4REC*(GAMMA-i.  )*UGMN*lGWN/2./AMN/AMN) 

QHEAT1*0.00298* (CV**0.333l ♦RHCMN*UGMN 
QHSAT»OEATl*(TGMN~TWAlL ) 

EGAS*-QHEAT  *PI  E*DI  A+AMGMMPCDPOT  *GZER0 
TAUW*0. 0005*RHGMN*UGMN*UGMN 
Fi=  Ur,MN*AMGMN  TAUW*PI£*CIA  4FD2A*AAEA 
Fl*Fl/770. 

RSl«AMGMN*<CV*TGMN»UGPN*UGMN/1556. ) ♦ < AREA*QEPLOT 
$*AMGMN/RHQMM*PSIMN/7?8. 

OSOT*AMGMN* (POOPOT*GZERC-( CV+RGAS/778* )*TGMN)  ♦ !TAl)W*UGMh/778.  - 


k. 


50 


200 


75 

300 


#QHEAT)*PIE*CIA4F024*UGFN*AR6A/778.-PSIMN*AREA*0EPI.DT/778. 
DS0T*DS0T/AREA/RHQMN/TGMN/ ( l.-EPLPN ) 

AREA*  ( ARE A2 4* ARE A22  ) / 2 • 

DELPOT=-ANN*UGMN*( AREA24-AREA22) / (X24-X22 ) /AREA 
0ELP0T*0ELPCT+778.*AMN*0SDT/RGAS 

WORK*(TAUW*PI£*OIA/AREA-AMN*RHOMN*DEPLDT+F0244UGMN*ANGMN/AREA)/RHO 

♦MN/U.-EPLMM 

WORK* WORK- ( AHN**3 ) *AHGPN/AREA/ < 1 .-EPLMN ) /PS IMN/GAMMA 

OELPOT*DELPCT-WORK 

P24*P22«-0ELP0T*(T24-T22) 

ANEW24*(P24-UG24)*(GAMFA-l.)/2. 

TG24* ( ANEW2 4*ANEW24 ) /G APPA/RGAS 
S24=S21*0SCT*(T24-T21)*778./RGAS/GAHHA 
EPL24*EPL21*D£PL0T*(T24-T21) 

IF ( EPL24.LE .0.00 )EPL24*0 .0 
WORK=  2.*GAPMA/(GAPMA— 1.  ) 

TG24* ( ANEW24*ANEW24 )/G*PPA/RGAS 
PSI24*P$I2l*l <AN£W24/A21M*W0RK) 

PSI24*P$I24*EXP(-GAHHA*(  S24-S21 ) ) 

RH024*PSI24/RGAS/TG24 

ABASE*BASE*(PSI24+PSI21) 

WORK* < ABASE/2. -R1)*GZER0/WT 
UGNE24»UG21*W0RK*( T24-T2 1 ) 

CHECK* ABS ( ( ANEW24-A24 I /A24 ) 

IF(CHECK.GT.TOL)  GO  TO  50 
IFUJG24.LE.0.C0)  GO  TO  50 
CH£CK*ABS(  ( UGNE24-UG24 ) / 1'624 ) 

IF (CHECK. GT .TOD  GO  TO  50 
WRITE(6*5000)QHEAT2*QHEAT1*QHEAT • EGAS* RSI 
XTRAV*(X24-0.125)U2. 

0EL=APGM6 I PODPOT*GZERC-CV«TGHN ) 

WRITE (6, 5000) L24, DEPICT, CSCT.AMGPKtDEL.XTRAV 

RETURN 

ITRY*ITRY«-l 

IFUTRY.LT.il)  GO  TO  2C0 
L0C*2 

CALL  01  AGO ( LGC*UGNE24*  ANEW24 ) 

UG24*UGNE24 
A24=ANEW24 
GO  TO  10 
X24*XLIP 

FQRPAT ( 1X.5E20. 5 ) 

T24*T22MX24~X22)/$L0PEP 

TMUZE*T24 

GO  TO  11 

ENC 

SUBROUTINE  RFORCE ( A, B ) 

01  PENSION  XX( 20 ) »RF ( 20  ) 

OATA  XX/O.C, C. 0333*0. 0416,0. 0499*0.066,0. 100*0. 1333*0. 

> 1666 f 0.200 *0.2333*0. 2666  *0.300*0 .333*0.500*0.6666* 1.00* 1.333* 1.50* 
$1,666*0,0/ 

DATA  RF/O.C, 520., 525., 520. ,420. ,255. ,215. ,195. *183., 17 

.,162. ,153. *14 5. *118. ,55. *52.* 1 8. * 7. 0*8. 0*0.0/ 


10 

15 


20 


A*A-0. 1249 
00  10  1*1,19 
IF(A.EQ.XX(  1) ) GO  TC  15 

IF(A.GT.XX(  I).ANO.A.LT.XX(Ul))  GC  TO  20 

CONTINUE 

8*PF( I ) 

6*8/2. 

A*AaO. 1249 
RETURN 

SLCPE*(RF(I«U-RF<n)/(XXU«U-XXU)) 


■jtcataMmi | 


ii 

a 


ii 


{ ! 

4 « 

ii 


li 


1 1 


30 


c»RFm-sioPE*xx<;i 

B**LGPE*A4C 

B-B/2. 

*•*♦0.1249 

RETURN 

END 

SUBROUTINE  BREACH 
CSUBROUT iNE  BREACH 

C THIS  SUBROUTINE  CALCULATES  THE  POINTS  LGCATEC  ON  THE  BREECH 
REAL  L5,L6,L3,L4,L21,L24,LlltL13 

C0PMQN/CGM2/X12,Tll ,Tl2,T13,UGll »UG12,UG13,A11,A12,A13,PSI11,PSI 12 
$»PSI13»TG11 ,TG12,TGL3,RHC11«RH012,RHG13,S11,S12,S13»EPL11*EP112,EP 
"L13.AREA11, AREA12*  AREA13  «L11,L13,LP12 
COPMGN/SAME/RGAS.PSIl, GAMMA, DZERC, PIE, BEE, PODPCT.CV, BASE, GZERC.WT, 
MOL, AN 

COMMON/EXIT /XLIM.TMUZE 

COMMON/ PARM /PAR  INI , PAR IN2, PARIN3 , PARBA1 , PARBA2 ,PARBA3,PARBRI , PARBR 
$2, PARBR3 

A1«SQRTIGAMMA*RGAS*530.) 

ITRY«1 

CONST— 0.58778E403 
TWALL-530. 

UG13-0.0 

AI3-A11 

RH0P«3.1085 

AREA* AREA 11 

Q12*2.*A12/ (GAMMA-1. ) 

PSU3*PSU1 

EPL13-EPLI1 

$L0PEGM-A13+UG12-Ai2>/2. 

T13—  X12/SLCPEQM12 
IF(T13.GE.TMUZE)  T13-TPU2E 
PSIMN* (PSI13+PSIll>/2. 

BR13*8EE*( ( PSIMN/144. ) **AN )/12 
0L13»BR13*(T13-T11) 

L13«L1I+0L13 

IF (L13.GT. 0.000625) L13«0. 000625 
CALL  DIMEILll iL 13, SURF 11 .SURE 13, VCLll, V0L13 ) 
IF(V0L11.LE.0.00I$URFIa*C.C0 
IF(VGL13.L£.0.00) SURF 1 3* C* CO 
£PLMN*(EPL134EPL11 ) /2. 

VOLMN«iVGLi 34VQL11 1/2. 

IF(VOLfctf.LE.O.OO)  GO  TC  31 
ETAMN-EPLMN/VOLMN 
GO  TO  32 
ETAHN*0.00 
CONTINUE 

AMGMN“ETAMN*8R13*RHGPM AREA 
AMGMN»AMGMN*( SURF1 3+SURF 1 1 1/2. 

DEPLDT— AMGMN/AREA/RHCF-EPLMN*UP12/X12 
FORMAT (1X«5E20. 5) 

TG13*A13*Al 3/GAMMA/RGAS 
TGMNS (TGi3*7Gil)/2. 

RHC13-PSI13/RGAS/TG13 
RHGMN» ( RhOl 3+RHO 11 ) /2. 

01  A* SORT (4. PARE A/PIE) 

FILM«3.24*EXP(CCN$T*T13)*0IA/2./Xl2 
QHEAT*FILM* (TG13-TNALL ) 

QHEAT«QHEAT*PARBR1 

EGAS— QHEAT*PIE*DIA+PCCPCT*AMGMN*GZERO 

R$l«AHGMN*<CV4TGMN4PSIMN/778./RHCMN)4P$INN*AREA*0EPLDT/?78. 
HORK“AREA*< i.-EPLMN )*RHQMN4TGHN 

DSOT»AMGMN* (PODPOT*GZERG-( CV4RGAS/776. )*TGMN ) -GHEAT*PIE*OI A - 


31 

32 


300 


HPSIMNMAREA4CEPLDT/778. 

DSDT«DSOT/hOAK 
IF(T13*GT.PARBR2)  GC  TC  200 
150  CONTINUE 

0£LQ0T»((Ali+Al3)/2.)*(778.40SDT/RGAS+DEPLDT/<n-EPLMN)) 
OELQOT«DELQCT*( (All/2* +A 13/2*)** 3 )* AHGNN/ ARE A/PS IMN/( 1*-EPLMN)/GAM 
AHA 

013«  Q1240ELQ0T*CT13-T12) 

P13-Q13 

ANEW 13* (GAMMA- l* )*Q13/2* 

TG13*ANEW13*ANEW13/GAMMA/RGAS 
S13*DSDT*( T13-T11 )*778*/RGAS/GAMHA+Sll 
PSI13*PSI11*( ( ANEW13/A11  >**(2.*GAMMA/(GAMhA-l.  ) ) ) 
PSU3*PSI13*EXP(-GAMMA*IS13-S11)  ) 

RHC13*PS I 13/RGAS/TG13 
EPLl3*EPLil*DEPL0T*(T13-Tll) 

IF(EPL13.LE.0*00JEPL13»0.0 
CHECK*ABS( ( ANEM13-A13) /A13 ) 

A13*ANEW13 

IF (T13*GE*TMUZE ) RETURN 
IF(CHECK*GT.TOU  GO  TO  180 
RETURN 

180  ITRY*ITRV+1 

IF ( ITRY.LT* 11 ) GO  TO  30 
L0C*3 

CALL  DIAGO(L3C,UG13,A13 ) 

RETURN 

200  WRITE (6,300 1 L 1 1 , L13 , VCL1 1, VGL13, SURF11 , SURF13, EPLMN, ETAMN, AM6MN, CE 
KPLOT  , QHEAT,EGAS«0S0T *AR£A 
GO  TC  150 
ENO 

SUBROUTINE  CIAGOILA, BUG, AA > 

C THIS  IS  A OIAGANOSTIC  SUBROUTINE 

REAL  L5,L6,L3,L4,L21,L24,L11,L13 

C0MM0N/C/X(99,99),TIME(9S»99),VEL(99*99)»SS(99,99)  , PRESS  (99,99)  , TE 
$MP(99»99)*DENS(99«99),VELP(99,99)  ,ENTR(99,99) , VOL FRA ( 99,991 ,ALI99, 
799) 

COMMQN/COMl /X3» X4»X6»T3,T4»T6,UG3 ,UG4,UG6, A3, A4, A6«PSI3, PS  1 4, PS  16, 
$TG3,TG4,TG6,RH03,RH04,RHC6 ,UP3*UP4,UP6,S3, S4, S6,£PL3, EPL4«EPL6,ARE 
>A3,ARCfl4,AREA6,L3*l4,L5,L6 

C0HH0N/CCM2/Xi2,Tll,T12,T13,UGll,LG12,UG13,All,A12,A13,PSUlfPSI12 
J,PSIl3,TGll,TGi2,TGl3,RHCll,RH012,RH0i3»Sll»S12,Si3,EPLil.fPH2,EP 
**L13,AREA11, ARE A 12, ARE A 13 ,Lil,L13,UP12 

C0HM']N/CCM3/X21,X22,X24,72l,T22,T24,UG2i,UG22,UG24,A2l,A22tA24,PSI 
S21,PSI22,PSI24,TG21,TG22  ,TG24 ,RHC21 ,RH022, RH024, S2 1, S22, S24,EPL2l , 
)EPL22»EPL24»AR£A2l»AREA22fAREA24,L21»L24,NSTART 

COHMON/SAHE/RGAS, PSI 1, GAMMA, OZERC, PIE, BEE, PODPOT,CV, BASE, GZERO,WT, 
>TQL,AN 

CQMMGN/CUT/MAXJ»MAXI 
GO  TO  (10, lit  12), LA 
10  WRITE46, 100 ) 

100  F0RMAT(1H1, 'EXCEEDED  ITERATION  LIHIT  IN  SUBROUTINE  INTER*) 


11  WR!TE<6,101> 

101  FORMAT ( 1H1 ♦ 'EXCEEDED  ITERATION  LIMIT  IN  SUBROUTINE  BASE*) 

GO  TC  15 

12  WR I TE( 6« 102  ? 

102  FORMAT ( 1H1 , ’EXCEEDEO  ITERATION  LIMIT  IN  SUBROUTINE  BREACH*) 

15  MR!TE(6,U0)BUG,AA 

110  FORMATUHO,  'LAST  CALCULATED  GAS  VELOCITY ', E12.3.//2X,  'LAST  CALCULA 
>TEO  SPEEO  Of  SOUND* ,£12*3) 

105  FORMAT (lX,I3,lX,I3,E12*3,2XfEl2*3,2X,E12«3«2X,E12*3,2X,E12*3,2X,Ei 
**2.3,2X,E12.3,2X,F6.3) 


maxi*maxi«i 
00  20  1*1. MAXI 
00  20  4*1,95 

IF (SSI I , J).EQ*0«0)  GO  TO  20 

WRITE(6,105)I«4,X(I.4),T1ME(1«J) * VEL  ( I . J ) , SS( I ,4) .PRESS! I , J) .TEMPI 
$1.4) »DENS! I »4),V0LFRA( 1. J) 

20  CONTINUE 
STOP 
■ ENO 

C SUBROUTINE  OUTPUT 
C THIS  SUBROUTINE  PRINTS  RESUTS 
SUBROUTINE  OUTPUT 

COMMON/C/ X( 99, 99), TIME (95, 99) .VEL (99,99) ,SS( 99,99) .PRESS (99.99 ) »TE 
*HP (99,99) , OEMS (99, 99 ) ,VELP (99,99 ) ,ENTR (99,99 ) , VOL FRA (99, 99 ) . AL( 99, 
?99 ) 

CONMON/SAME/AGAS, PSIl, GAMMA, OZERC, PIE, BEE, POOPCT.CV, BASE, GZERO.WT, 
MOL  ♦ AN 

C0MM0N/0UT/MAX4,MAXI 
MRITE(6,100 ) 

100  FORMAT ( 1H1 , •**•**♦  SOLUTION  TO  GAS  FLOW  IN  TUBE  *********•) 

WRITE (6, 104) 

104  FORMAT (lhO, 'FLOW  PARAMETERS' ,//3X,' l * »3X, 4 4' ,3X, 'POSITION* , 8X, 'TIM 
>£', 9X, ' VELOCITY', 3X, 'SPEED  OF  SOUND* ,3X, 'PRESS' ,3X, 'TEMP* ,3X, 'DENS 
$ITY»,3X, 'FRACTION') 

105  FORMAT (IX, 13, IX ,13, E12«3,2X,£12«  3,2X,E12«3»2X«E12*3,2X,E12»3,2X,E1 
«2.3,2X,E12.3,2X,F6.3) 

MAXI  I *MAXI ♦ 1 
00  20  I*1,MAXII 
00  20  J*l ,95 

IF(SS(I,4).EQ.0.0)  GO  TO  20 
PR PS I “PRESS ( I , 4 )/l44« 

WR ITE( 6,105  )I ,4,X(I,J),TIME(I,J),VEL(I,4) , SS ( I , J ) , PRPSI ,TEMP ( I, J ) , 
SDENS ( I , J ) , VCLFRA( I , J ) 

20  CONTINUE 
WRIT£(6,110) 

00  30  I*1,MAXII 
00  30  J*l,95 
I F ( I .NE.JJ  GO  TO  30 
XTRAV*(X(I, J)-0.125)*12. 

PRPSI *PRESS ( I,J)/144# 

WRITE(6,1U)I,J,X(I,J),TIME(I,J)  , XTRAV,  PRPS I , VEL  ( I , J ) 

30  CONTINUE 

110  FORMAT (1H1.20X, 'PROJECTILE  DATA* ,//3X, ' I ', 3X, • J' , 3X, ' POSITION' , 8X, 
• 'TIME', 9X, 'TRAVEL', 8X,  'PRESSURE' ,8X, 'VELOCITY' ) 

111  FQRMATdX,  1 3,1X,  13,  E12«3,2X,E12«  3,2X,F12«2,?X,  F12«0,2X,F12*0) 
RETURN 

ENO 

SUBROUTINE  PROP(A) 

C SUBROUTINE  FCR  CALCULATING  PROPERTIES  AS  A FUNCTION  OF 
C FRACTION  OF  DETERREO  ANO  NCNCETERREC 

COMMON/SAME /RGAS, PSI 1, GAMM A, OZERC, PIE, BEE, POOPOT.LV, BASE .GZERO.WT, 
>TOL, AN 

C0MM0N/PCW0ER/CV1.CV2 , RGAS1,RGAS2 ,GAM1,GAM2, OLAY 
A*A* 12* 

AMN0*2.0175E-0?4(A-0LAY*12.)/(0.0CT5-0LAV*12.) 

AMT*6.815E-08+AMN0 

AMFNO*AMNO/ AMT 

CV*CV1*( 1«*AMFN0) ♦CV2*AMFNC 

RGA$*RGASl«(l.-AMFND)tRGAS2«AMFNC 

G~NMA*l.*RGAS/CV/778. 

A*A/12, 

RETURN 

ENO 
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C SUBROUTINE  START 

REAL  L5,L6,L3,L4,L2l,L24,LU,L13 

CONMON/CCM2/X12»Tll,T12,T13,UGll,LG12.UG13,All,Al2,A13,PSIli*PSll2 
A«PSI13»TG11*TG12»TG13*  RHCl 1 ,RH0l 2 »RH0I3» Sll»Sl2*Sl3*EPl 11,6PL12.EP 
$Ll3»AREAll»  AREAl2»AREA13iLlltLl3fUPI2 

COMMON/SAME/RGAS»PSI1 » GAMMA, DZERC, PIE* BEE,  PODPOT,CV, BASE, GZERO,WT, 
ITOL* AN 

COMMON/ STAR /K ST ART 

AK*KSTART 

ITRVi 

UG13-0.0 

AI3*A12 

DELX*-0. 125/AK 

P$113*PS112 

RHGP*3.1085 

X13*X12+DELX 

AREA12*TAREA(X12) 

AREA13*TAREA<  X13) 

AREA«(AREA13*AREA12)/2. 

Q12»2.*A12/ (GAMMA-l. ) 

EPL13*EPL12 

30  SLOPEQ*-! A13+A12J/2. 

PSIMN*(PSI 12+PSI 13) /2. 

T13*T12+ ( X13-X12I /SLOPED 
BR13»BEE*UPSIMN/144.)A*AN)/12> 

0L13*BR13*(T13-T12> 

L13«L11*0L13 

CALL  DIME  (LU  ,L13«  SURF12  ,SURF13*  VCL12,VQL13) 

EPLMN*(EPLl3*EPLl2)/2. 

V0LMN*lVQL13+V0L12)/2. 

ETAMN*EPLMN/VOLMN 

AMGMN»ETAMN*BR13*RH0P*AREA*<SURF13ASURF12 ) /2* 
DEPLDT*-AMGMN/AREA/RHCF 
TG13*A13*A1 3/GAMMA/RGAS 
TGMN*  <TG13*TG12)/2. 

RHC13*PSI 1 3/RGAS/TG13 
RHCMN* ( RH01 3+RH012 ) /2. 

EGAS*POCPGT*AMGMN*GZERC 

RSl*AMGMN*tCV*TGMN+PSIPN/778./RHCPN)^PSIMN*AREA*DEPLDT/778. 
WORK*AREA*( l.-EPLMN)*RHCMN*TGMN 
DSOT*<  EGAS-RSI ) /WORK 

OELQOT*(  ( A12+A13)/2« )*  (778.*DSDT  /RGAS+DEPLCT / ( l.-EPLMN) ) 
DELQDT*DELQCT+< 1 ( A12+A 13  )/2. )**3 ) 4AMGMN/ AREA/PS 1MN/ 1 l.-EPLMN) /GAMM 
#A 

Q13*G12*CELC0T*<T13-T12) 

P13*013 

ANEW13*(GAMMA-1. )*Q13/2. 

TG 13* ANEW 1 3*ANEHi3 /GAMMA /RG AS 
S13*S12*CS0T*778.M(T13-T12 )/R GAS /GAMMA 
PSri3*PSI12A( (ANEW 13/ A 12 )♦* (2.*G AMMA/ 1 GANMA-1* ) ) ) 
PSU3»PSU3*EXP(-GAMMA*(S13-S12)  ) 

RH013*PS1 13/RGAS/TG13 
EPU3*EPL12*0EPL0T*<T13-T12) 

CH£CK*ABSUANEWl3-A13)/A13) 

Al  3*ANEW13 

IF1CHECK.GT.T0U  GO  TO  180 

X12*X13 

RETURN 

180  I TRY* l TR V>  1 

lFUTRY.LT.il)  GO  TO  3C 
L0C*3 

CALL  01 AGO! L€C, UG13 , Al 3 ) 


RETURN 

ENO 

FUNCTION  TAREAULOC) 

C THIS  IS  A SUB  PROGRAM  TO  CALCULATE  TEE  BORE  CROSS  SECTIONAL  AREA 
DIMENSION  BCRE( II ) » POS (11) 

DATA  BORE /0.024900, 0.0 2491, 0.0 2771 ,0.02746, 0.02514, 0.01864, 0.018 

>64, 0.01824,0.01824, 0.01824,0.01824/ 

DATA  POS/-.10,0.0, 0.0279?, 0.1032,0. 1078,0. 1135, 0.1300, 1.542. 

$3.0,5.0,15.0/ 

DO  10  1*1,11 

IFiPGSin.EC.ALCC)  GO  TO  5 
IF(POSm.GT.ALOC)  GO  TO  7 
GO  TO  10 

5 TAREA*3.141593*B0R£ ( I >*BCRE( I )*0 .25 
GO  TO  15 

7 AM«<BOREU-l)-BCREU)l/<POS<I-i)-POSUn 
B*BORE ( I ) 

00*B+AM*ULCC-PGS(in 
T AREA* 3. 141593*00*00/4. 

GO  TO  15 
10  CONTINUE 
15  RETURN 
END 

SUBROUTINE  CIME1A,B,SUPA,SUR8*VOLA,VOLB> 

CSUBRCUTINE  DIME 

C THIS  IS  A SUBROUTINE  TO  CALCULATE  PROPELLANT  SURFACE 
C AREA  AND  VOLUME 

DATA  PIE, ABAR,R 1/3. 141593, 0.00095*0.0006249/ 

WORK=RI-A 

WORKl*RBAR-RI 

VOLA»2.*PIE*WORKl*WORKl*WORMPIE*PIE*WORK*WORK*<WCRKl4 

>4.*W0RK/3./PIE) 

SURA«2.*PIE*(PIE*U0RK*(WCRK142.*NCRK/PIE)4U0RK1*M0RKI) 

W0RK*RI-8 

V0LB*2.*PIE*WQRK1*WGRK 1* WORK*PIE*PI E*WORK*WQRK*( W0RK1* 

>4.*W0RK/3./PIE) 

SURB»2.*PIE*(PIE*W0RK* (WCRK1+2.*WCRK/P IE  )+MORKl*WCRKl ) 

RETURN 

ENO 

SUBROUTINE  INTER 
C SUBROUTINE  INTER 

C THIS  SUBROUTINE  CALCULATES  INTERIOR  NETWORK  POINTS 

C 

C 

REAL  L5,L6,L3,L4,L21,L24,LU,L13 

C0MM0N/CGM1 /X3 , X4,X6»T3»T4»T6,UG3  ?UG4,UG6, A3,A4, A6,PSI 3, PSI  4,PSI6, 
$TG3,TG4,TG6 ,RH03,RH04, RHC6,UP3,UP4,UP6, S3, S4, S6,EPL3,EPL4,EPL6, ARE 
>A3,AREA4,AREA6,L3,L4,LS,L6 

COMMON/SAME /RGAS, PSI 1* GAMMA* 0ZERC*P IE, BEE, POCPOT,CVf BASE tGZEROtWT, 

>TOL , AN 

COMMON/EX  IT /XL  I M,T MU ZE 
C0MM0N/8HEAT/FACT,TWALL,FR!CT 

COKMON/PARM/PARIN1,PARIN2,PARIN3,PARBAI,PARBA2,PARBA3,PARBR1,PARBR 
$2 , PARBR3 

DATA  AMU, REC,RHOP/ 1.0445 £-06*0.9, 3. 1085/ 

FUG*1.0 

U»SQRT(GAMMA*RGAS*530.) 

1TRY*1 

UG6*IUG3*UG4)/2. 

A6* ( A3+A4 ) /2. 

TG6«ITG34TG4)/2. 

RHC6* I RHQ3*RK04 ) /2 • 

PS  1 6* I PSI 3* PS 14 1/2. 
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P3*2.4A3/CGAMHA-1.)  +UG3 
04*2.4A4/(GAMMA-i. )-UG4 
6PL6*(EPL3+EPL4 )/2. 

AREA* (ARE A3* ARE A4) /2. 

S6*($3^S4)/2. 

C CALCULATION  OF  THE  LOCATION  CF  THE  NEK  POINT  ON  X-T  DIAGRAM 

30  SL0PEPc(UG3*A3*UG6+A6)/2. 

$L0PEQ«(UG4-A4*UG6-A6)/2. 

T6*( X4-X3+T34SL0PEP-T44SL0PEQI/C SLOPEP-SLOPEQ ) 
IF(T6.GE.TML2E)  GO  TO  25C 
X6*X3+SLCP£P4<T6-T3) 

31  AREA6*TAREA (X6) 

K0RK*X4-X3 

IF(T4.EC.T3)  GO  TO  32 
T7»X6-X3^T3*W0RK/(T4-T3)-T6*UG6 
T7*T7/ ( WCRK/( T4-T3 l-UGt ) 

GO  TO  33 

32  T7*T3 

33  CONTINUE 
X7*X6-UG6#(T6-T7I 

IF (T4.EQ.T3  ) GO  TO  34 
T5*X6-X3+T3*K0RK/ C T4-T3 ) -T6*UP6 
T5*T5/ (WORK/I T4-T3 J-UP6 ) 

GO  TO  35 

34  T5*T3 

35  CONTINUE 

X5*X6— UP6*( T6-T5) 

IF(X5.LT.0.C>  X5*X3 
WORK l* ( X7-X3) /WORK 
UG7*UG3+(UG4-UG3)*W0RK1 
A7*A3+ ( A4-A3) •WCRKl 
S7*S3+ ( S4-S31 *WGRK1 
TG7»A7*A7/<GAMMA*RGAS) 

W0RK2*2.'*GAPMA/ (GAMMA-1.  ) 

PSI7«PSI3*( (A7/A3)*4M0RK2J 
PS17»PSI7’*E>P(-CAMMA*( S7-S3I ) 

RH07*PSI7/RGAS/TG7 
WORK!* (X5-X3) /WORK 
UG5*UG3*-  ( UG4-UG3  >*WORK  1 
UP5*UP3+(UP4-UP3)*H0RK1 
A5*A3+( A4-A3) ♦WGRK1 
$5*S3+ ( S4-S3 ) *W0RK1 
TG5*A5*A5/ ( GAMMA*RGA$ ) 

P$I5*PSI3*( (A5/A3) **H0RK2I 
PS15*?SI5*EXP(-GAMMA4( S5-S3I ) 

RHC5*P$I5/RGA$/TG5 

EPL5*EPL3*(EPL4-EPL3)*KCKK1 

IF(EPL5.LE.G.00)EPL5*0.0 

L5*L3ML4-L3)4MGRK1 

UMN*(UG6tUG7)/2. 

UGMN*UMN 
AMN* ( A6* A7 ) /?. 

TGMN«(TG6+TG7)/2. 

RHOMN* ( RH06>RH07 ) /2 • 

PSIMN* (PSI6-*PSI7)/2. 

UP56* (UP5+UP6 ) /2* 

RELV5*ABS(UG5-UP5) 

RELV6*A8$(UG6-UP6) 

RELV56*(RELV5*-RELV6)/2* 

PEE56*(PSI5*PSI6)/2. 

BR56*BEEA( (PEE56/144. JA*4N)/12. 

IF  RELV56.GT.500.I  BR5e-eR56*Il.*C*Q004*<RELV56-500. )) 


DL56*BR56*I T6-T5 ) 

L6*L5+DL56 

IF IL5*GT .0.CGQ625) 15*0 .CC0625 
IF ( L6*GT«0*C00625 ) L6=0  .0C0625 
CALL  DIME<L5»L6»SURF5»SURFfc,VCL5tV0L6) 

IF (VGL5.LE. 0.00 ) SURF 5*C«C0 
IF ( V0L6*LE«G*00 ) SURF6*C  »C0 
EPL56«(EPL6<EPL5)/2. 

0BAR*DZ6R0- L5-L6 
IF (OBAR.LT. C. 000652 )DBAR*0.00 
REY56«RELV56^DBAR/AMU 
REY56=REY56*( RHQ5+RH06 1/2* 

IF (REY56«LE *0*000)  GO  TO  98 
C0RAG*28*/( REY56**(0*85) WC*48 
GO  TO  99 

98  C0RAG«0.0 
GO  TO  99 

C CALCULATION  CF  VELOCITY  ANC  PRESSURE  GRADIENTS 

99  IF(T3.GT.T4)  GO  TO  100 
T8*T4 

X8=(UG3+A3)*<T8-T3>+X3 

PSI8*PSI3'M  PSI6-PSI3)*IX8-X3)/("6-X3 ) 

UP8*UP3+  ( UP6**UP3  )*(  X8-X3  )✓  (X6-X3  ) 

0P0X«0.00 

0UP0X*(UP4-UP8)/(X4-X8> 

GO  TO  125 
100  T8=T3 

X8* ( UG4-A4 )*(T8-T4)+X4 
PSI8*PSI4*IPSI6-PSI4)*(X4-X8)/(X4-X6) 

UP8*UP4* ( UP6-UP4 )*  < X4-X8  )/ ( X4-X6  ) 

0P0X«0.00 

0UPDX*(UP8-LP3)/<X8-X3) 

125  CONTINUE 

C ORAG  FORCE  ON  PROPELLANT  PARTICLE 

F0R56=CDRAG*{ RHC5+RH06 )*RELV56*RELV56 
'F0R56*FCR56«»DBAR*08AR/16. 
F0R56=F0R56-DP0X*<DBAR**3)/6. 

F0R56=FCR56*P1E 

IF(VCL5.LE.C.00.AND.VCL6 .LE.0«00 > GO  TO  1125 
ETA56=2.*EPL56/(V0L5»VCL6> 

GO  TO  1126 

1125  ETA56*0.00 

1126  CONTINUE 
FD56«ETA56*F0R56 

C MASS  OF  GAS  GENERATED  GOING  FROM  5 TC  6 
AMG56=ETA56*8R56*RHCP**REA 
AMG56*AMG56*( SURF5  + SURF6  »/2. 

AMG56»AMG56*FUG 

C CHANGE  IN  PARTICLE  VELGCITY  ANC  MASS  FRACTION 
IF(EPL56*LE*0»  >000)  GC  TC  126 
DUPOT  »UP56*AMG56/RHCP/ AREA/EPL56 
OUPD  r*OUPOT ♦F056/RHCP/EPL56 
GO  TO  127 

126  0UP0T«0.0O 

127  CONTINUE 

06PLDT»-EPL56*0UPDT/UP56 
DEPL0T*0EPLCT-AMG56/RHCP/AREA 
01  A* SORT (4«*ARE A/PIE) 

FI LM=FRICT*GAMMA*RGAS*RHCMN*UGMN/2*/ (GAMMA-i# )/778* 
TAW*TGMN*( 1 • ♦REC*<  GAMMA- 1«  )*UGMN*UGMN/2» /AMN/AMN ) 
CHEAT«FILM*ITAM-T*ALL) 

QHEAT=CHEAT*FACT 

EGAS=-0HEAT*PIE*0IA-F056*AREA4RELV56/778. 


!i  c 


bWMd«cuMO'rMrt)}0*KUUKUI  *6/tK(j 

TAUW*0*5*FR  ICT*RHOMN*UGPN*UGMN 

FI*  UGNN*AMG56*TAUW*PI£6CIA*F056*AREA 

Fl*Fl/778. 

RS1*AMG56*(  CV*TGMN*UGPMUGPN/1556*)*PSIMN*AREA*DEPL0T/778* 
R$1«RS1*PSIHN*AMG567RHCHN7778. 

WORK-AREA*( 1.-EPL56 )*RFOPN*TGMN 

0SDT*AMG56* (PQOPOT*GZERO~( CV+RGAS7778* )*TGMN ) ♦FC56*AREA*(UGMN-RELV 
156)7778*  ♦ ( TAUW*UGMN/778 .-CHEAT ) *PIE*Oi A -PS IHN*AREA*0EPLDT7778* 
DSOr-OSDT/WCRK 
IF (T6*GE*PARINl ) GO  TO  201 
128  CONTINUE 

CALCULATE  THE  CHANGE  IN  P ANC  C 
UG36* (UG3+UG6 ) 72* 

RHQ36* ( RH03 ♦RH06 ) 72* 

PS!36*{P$I3+?SI6)/2* 

AREA*  ( AREA3+AREA6 ) 72* 

A36* ( A3+A6) 72* 

EPL36*( EPL3*EPL6)72* 

TRY»A36*A36*778**<  S6-S3 ) /( X6-X3 ) 7CAMHA7RGAS 
CELPDT«-A36*UG36*(AREA3~4REA6)7m-X6)7AREA 
0ELPDT*0ELPCT*A36*778.*0$0T7RGAS 
WORK*TAUW*P  IE*OI A/AREA7 < 1.-EPL36  17RH036 
W0RK*WCRK+FC567RH0367( 1.-EPL36) 

DELPCT*DELPCT-WORK 

OELPOT*OELPCT ♦A36*0EPLCT /< l.-EPL  36 )+( ( A36**3 )/PSI 36/GAMMA-UG36/RH0 
>36 ) *AMG567AR£ A7 ( 1.-EPL36 ) 

DELPCT*OELPCT*TRY 
PS  146*  (PSI4*PSI6)72* 

EPL46* ( EPL4+EPL6 ) 72* 

A46«(A4*A6) 72. 

UG46* ( UG4+UG6 ) 72* 

RH046= ( RH04+RHQ6 ) /2. 

AREA* ( AREA4+AREA6 >72* 

TRY»A46*A46*778.*I S6-S4 ) /( X6-X4 ) 7GANHA7RGAS 
0ELQDT»-A46*UG46*( AREA4-AREA6)/! X4-X6)/ARE/ 
0ELQ0T*DELQ0T+A46*778.*CSCT7RGAS 
WORK*! TAUW*PIE*OIA/AREA+F056 >7(1 *-EPL46)7RH046 
DELQCT«DELQCT*W0RKtA46*0EPLCT71 1 .-EPL46) 

OELGDT*OELQCT*! 1 A46**3 )7PS I46/GAPPAHIG467RH046 )*AHG567AREA7 U.-EPL 
146) 

DELCDT*DElCCT-TRY 

P6*P3+0ELP0T*!T6-T3) 

Q6*G4*0EL001*!T6-T4) 

ANEW6*IGAWMA-1. )*<P6*Q6)74. 

UGNEW6*(P6-C6)/2. 

S6r$7*0S0T* (T6-T7 )* 778  */RG AS/GAM PA 
EPL6*EPL5+0£PIDT*CT6-T5> 

IF ( EPL6*LE* 0*00 )EPL 6*0*0 
UPN6»UP5tDUPDT* ( T6-T5 ) 

IF(EPL6*LE  *0*00 )UPN6*LGNEb6 
TG6*ANEW6*AhEW67GAMMA/RGAS 
P$I6- P5I ?*t  UNEK6/A-;  }**WCRX2) 

PSI6*P$I6*EXP!-GAPP.A*I$6-S7) ) 

RH06*PSI6/RGAS/TG6 
I F ( T6*GE*TMtZE ) GO  TO  260 
5003  F0RMAT<lX*5Eli.3) 

CHECK*A8S( ( ANEW6-A6 )/A6 ) 

I F (CHECK *GT *TOL ) GO  TO  180 
CHECK*A9S( ( LGNCW6-UG6 ) 7UG6 ) 

IF(UG6*I T.3C0.)  CHECK-C.C0C01 
TOLl*lO.*TCL 

IF(CHECK*GT .TOLL)  GO  TC  180 


UQLD4*UG4 
A0LD4*A4 
X0LD4*X4 
S0L04-S4 
PSIQL4*PSI4 
UP0L04*UP4 
XCL04*X4 
RETURN 
ITRYMTRYU 

IFUTRY.LT. 21)  GO  TC  2CO 
LOC«l 

CALL  OIAGG( LCC»UGNEW6» ANEW6 ) 

UG6*(UGN£W6+UG6)/2. 

A6« ( ANEW6+ A 6 ) /2  « 

UP6» ( UPN6+UP6 ) /2. 

GO  TO  30 

WRITEI6* 5004 ) V0L5»V0L6 *SLRF5»SURF6*ETA56t AMG56*EPL56»DUP0T*O£PL0T 
WRITE(6,5004)L5,L6,X5,X6»BR56 
WR I TE  ( 6 * 5004 ) P6  * Q6  * A6  » LG  6 « 0P6 
OVOL*2.4( < 3.*V0L5/ (4.4F I E ) )**0.33) 

WRITE (6*5004 )F056*TAUW*EGAS*RSl* F 1* UGMN*  F0R56* UP3*UP4* UP5*0UP0X * UP 
| M 156 *UG5 *OBAR  *0V0L»UPN6 

| 1 5004  FORMAT (1X*5E20. 5) 

f GO  TC  128 

M 2f 0 T6-TNUZE 
If:  X6*X3+SLCPEP*(T6-T3) 

»AQ« ( UG3*U0lC4-A3*A0LD4 )/2. 

AP*( UOL04+UG4+A4+AQLO4 ) / 2 • 

* I TI«(X4-X6*AC*T6-AP*T4)/(AQ-AP) 

I XI »X4+AP4( T I-T4 I 

[ \ Uo4»(UG4-U0t04)*(Xl-X4)/<XCLD4-X4)*UG4 

A4* ( A4-A0L C4 ) * ( X I-X4 1 / ( X0L04-X4 ) +A4 
; 1 TG4>A4*A4/GAHHA/RGAS 

PS14»(PSI4-PSI0L4)4(XI-X4)/(X0L04-X4)+PSI4 
f I RHQ4-PSI4/RGAS/TG4 

? \ UP4«(UP4-UPCL04)*(XI-X4)/(XCL04-X4)4UP4 

f } S4* ( 54-SGL04)6(XI-X4) / (XCLC4-X4 ) +S4 

;•  | X4-XI 

t ( T4*T I 

GO  TC  31 
260  UG6«UGNEW6 
\ , A6-ANEW6 

f ? UP6-UPN6 

RETURN 
ENO 


1 t ]\  JOB  MCASSEY  6044999919 
| /j j exec  watftv 


/PROGRAM  MCASSEY 

C PROGRAM  CALCULATES  THE  PROPAGATION  OF  A FLAMEFRONT  THRU  A 
C POROUS  PROPELLANT  BED 

C F-PROPELLANT  IMPETUS  (FT-L8S/L3) 

C C * CHARGE  WT  (LBS) 

C RHOP*  PROPELLANT  TENSITY  (LBS/IN**3) 

C ALFA«FLAMEFRONT  VELOCITY  PRESSURE  COEFFICIENT  ( IN/SEC-P$I**BETAI 

C BE'’ a*  FLAME  FRONT  VELOCITY  PRESSURE  EXPONENT  (0IM6NI 

C B*  PROPELL  ANT  BURNING  RATE  COEFFICIENT  HN/{  SEC-PSI**EN)I 
C CN* PROPELL ANT  BT  - NG  RATE  EXPONENT  (DIMEN) 

C RO=EQUIV.  MEAN  K JUS  OF  PROP.  GRAIN  (IN) 

C R 1=  1/2  WEB  ' TN) 

C XLO=LENGTH  OF  COMPACTION  REGION  (IN) 

C ETA-PROPELLANT  COVOLUME  (IN**3/L8) 

01  MENS  I ON  POO  ) ,PBAR(30)  * VF <30 ) ,DL ( 301 ,XL(30,30>,RL(:0,33)  »V(  30*30 
$)  »XMG(30»30)»ZL(30)  ,C».  (30)*E(30»30) 

C SUBSCRIPT  J DENOTES  THE  NUMBER  OF  TIME  INTERVALS  (PT)  DURING  WHICH 
C THE  FLAMFFRONT  PENETRATES  THF  PROP  ELL  ANT , J- 1 CORRF  SPONDE  TQ  IGNITIO 
DATA  RHOP ,°0,R I* XLO, ETA/O. 0560, 0.0115, 0.0075 ,0.23, 32. 2/ 

2 READ(5,5000)F, ALFA, BETA, B,fcN,. -ACT 
REAC( 5, 5000 )RGAS, GAMMA, CV,TZiRO 
5000  FORMAT ( 2E20. 5 ) 

WR ITE ( 6 » 5000 ) F , S LFA, BETA ,B,EN, FACT 
W R I T F ( 6,5000)RGAS,GAMMA,CV,TZERQ 
VEL*0.0 

AZERO«SQRT(GAMMA*RGAS*TZERO) 

C-4.014E— 03  * 

I COUNT-O 
00  10  1*1,30 
00  10  J = 1 * 30 
F ( I , JJ-0.0 
10  XL( I , J ) *0.0 
PIE-3. 1415 
T =0.0001 
XTG-0.0 
PRsRO-R  I 

V0I*2.*PIE*RP*»R*RH-P!e*PIE*RI*RI*(Rfi*1.333*RI/Pie) 

C XN* THE  NUMB EP  or  GRAINS  PER  UNIT  VOLUME 
C A-CROSS  SECTIONAL  AREA  OF  CARTRIDGE 
C 0 = 01  AMFTfr  OF  CARTRIDGE  C.ASF  0.335 
A=RIE*n. 335*0. 335/4. 

C VC* VOL(JMc  OCCUPIED  BY  COMPRESSED  PROPELLANT 
Vn*A*( 1.5-XLO) 

XN-C/ :R«OP*VOI*vC  ) 

*»R  I TE  ( 6 * 1234 ) 

P ( 1 ) *2  750. 
iv'i  1000  J*2 , 29 
P(J)«P( J-i) 

OT=0.?5E-04 

TtT*-DT 

H0\  I COUNT* If OUNT ♦ 1 

C SUBSCRIPT  j DENOTES  THE  NUMBER  OF  TIME  INTFRVAL(DT)  DURING  WHICH 

; T'<H  SLAMEF«GNT  PENETRATES  THE  PROPELLANT  3E0IIGNIT  T I M E* 0. 0001  S ) 


> A L C H A Tc  ThE  DISTANCE  (O(J)J  THE  FLAMFFRONT  HAS  TRAVELED 

PrtAR/J)-(P(J-i)*P{j 

VR  ( J )*AL?A4PBAR  ,•  j )**BE  TA 

OK  J)*Vp(  J)*DT 

ZL ( 1 ) * XLO 

Cl  ( 1 ) * X L n 


.*  38 


IF<ZL<  JJ.GE.1.5)  ZUJ)«1.5  <7 

IFCZLI J).GE.l.5)DL< J)*ZLtJ)-ZL<J-l ) 

IF(ZL( J).GE.1.5)T*T-0T 
IF(ZL( J).GE.1.5)0T«DL( J)/VF(J) 

IF(ZLCJ).Ge.l.5)T»T+0T 

CALCULATE  THE  DISTANCEULC  1,JI)  BURNED  NORMAL  TO  PROP.  SURFACE(IN) 

CELL  I DURING  TIME  OT,I*l  IS  INITIAL  COMPACTION  REGIONINO  PROP. I 

PXTG-0.0 

VOL  S*0. 0 

no  500  I«2» J 

0XL*B*P8AR(  J I **EN*DT 

DXL«DXL*FACT 

XL(I»J)*XL(I»J-1I+DXL 

RL(I» J)*PI-XL( It J) 

V(I,J)*2.*PIE*(RR**2*RLt  I , J ) HP  IE*  PI  E*RL(  I , J )**2*( RRH  .333/P I E*RL  ( 
#I»JJ) 

XMG( I y J ) IS  THE  AMOUNT  OF  GAS  GENERATED  FROM  CELL  I DURING  TI ME  T 
XMG ( I * Jl *DL ( I J *XN*A* ( VOl-V  < {» J ) >*R  HOP 
XTG=THE  TOTAL  GAS  GENERATED  IN  TIME  T 
nXMG=XMG(  i, J) 

EU,JJ*(XN*A*OLm*VOI-XMGU»J>/RHOPI/A/PUn 
VOL$«VOLS+A*DL ( I ) *6  U , J ) 

PXTG-PXTG+DXMG 
500  CONTINUE 

XNUM* 1 2*  *PXTG*F+P ( 1 ) *A*ZL ( 1 ) 
n£NOM»A*ZL ( J ) -VOL  S-PXTG*ET  A 
PPB  AR*  XNUM/ DENOM 
I MN*J 

I F ( ABSi I PPBAR-PC J ) ) /rPBAR) -0.05)55  » 55 » 800 
800  TF< IC0UNT.GT.25)  GO  TO  1001 
P( J)«PPBAR 

GO  TO  801 
55  CONTINUE 
XTG-PXTG 
P( J)*PPBAR 
I COUNT  -0 

nn  900  i =2 » j 

W°ITE(6»2468)T ,P( J) , ZL ( J I , XTG, P ( I , J) , XL( I. JJ,I,J 
TIME-T 

900  CONT4NUE 

I F ( ZL ( J ) • GE .1 . 5 ) GO  TO  1001 
K*J 

CL ( J >*ZL< J > 

3234  FORMAT (lHlt9Xt* TIME* ,15X,*HP<J) *,15X, • ZL { J > » ♦ 1 5X, « XTG* , 1 5X. *E (I t J I 

*'  ) 

2468  FORMAT (6( 4X  » £ i 5 . 6 ) » 2 1 3 i 

1000  continuf 

00  1002  J *?  f 9 
nn  1002  1*2,9 

wRIT6(6,2468)XLU  , J ) ,DXL , V II  , J ) , XMG<  I » J ) , RL  (I  , J ) »DL(I>  » I t J 

i002  continue 

GO  TO  1005 
1C01  CONTINUE 
J ~ I MN 

p(  JJ-*0(  J)*144. 

nE\c>=P(  I)  /PC, AS/ TZ  PRO 

“ N r R * {RGAS/ 77d.*CV)*5LOG(  T ZERO/530.  l-P.GAS/ 778. *ALOG<  P<  J)  /20  75.  ) 
c NTRsPNTR*778. /RGAS/GAMMA 
CO  1004  I-l.TMN 
)!GT*ZL( I 1/12. 

XuC I » J ) *XL( I tJJ/12. 


92  TIH«TM0.125-OI$T)/AZERO  ** 

93  WRITE (7,4445) I, J 

94  WR1TE(7,4444)DI$T,T!M,E(I , J ) , AZERO  , VFl  ,Vi=l,P{  J » ,DENS,E  NTR,  XL  (I  , J) 
S,TZERO 

95  WR ITE( 6, 4444)0 IST,TIM,E(I»J) » AZERQ,  VFL *VEL,P(J) »D6N$»  ENTR , XL ( ! * J) 
$»  T ZERO 

96  WRITE! 7,4445 ) 1 , J 

97  1004  CQNTINUP 

98  4445  F0RMAT(5X, •***♦****•, 5X, 14, 5X , 1 4, 5X, •*********' 5X ) 

99  4444  F0RMATUE20.5) 

100  1005  CONTINUE 

101  GO  TO  2 

102  END 


/GO 

0*  32060E  06 
0.48000E  00 
0.70000E  00 
0.25643E  04 
0.11470E  02 


0.31500E  02 
9.26300E-02 
0.10000E  01 
0.12506E  01 
0.46150E  04 


